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Abstract 

Sensor  noise  is  an  unavoidable  fact  of  life  when  it  comes  to  measurements  on  phys¬ 
ical  systems,  as  is  the  case  in  feedback  control.  Therefore,  it  must  be  properly  addressed 
during  dynamic  system  identification.  In  this  work,  a  novel  approach  is  developed  toward 
the  treatment  of  measurement  noise  in  dynamical  systems.  This  approach  hinges  on  proper 
stochastic  modeling,  and  it  can  be  adapted  easily  to  many  different  scenarios,  where  it 
yields  consistently  good  parameter  estimates.  The  Generalized  Minimum  Variance  algo¬ 
rithm  developed  and  used  in  this  work  is  based  on  the  theory  behind  the  minimum  variance 
identification  process,  and  the  estimate  produced  is  a  fixed  point  of  a  mapping  based  on  the 
minimum  variance  solution.  Additionally,  the  algorithm  yields  an  accurate  prediction  of  the 
estimation  error.  This  algorithm  is  applied  to  many  different  noise  models  associated  with 
three  basic  identification  problems.  First,  continuous-time  systems  are  identified  using  fre¬ 
quency  domain  measurements.  Next,  a  discrete-time  plant  is  identified  using  discrete-time 
measurements.  Finally,  the  physical  parameters  of  a  continuous-time  plant  are  identified 
using  sampled  measurements  of  the  continuous-time  input  and  output.  Validation  of  the 
estimates  is  performed  correctly,  and  the  results  are  compared  with  other,  more  common, 
identification  algorithms.  The  Generalized  Minimum  Variance  results  are  generally  better 
than  those  of  the  other  methods. 
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COUNTERING  THE  EFEECTS  OF  MEASUREMENT  NOISE  DURING  THE 
IDENTIFICATION  OF  DYNAMICAL  SYSTEMS 

I.  Introduction 

System  identification  entails  an  empirical,  data  driven  approach  to  modeling  dynam¬ 
ical  systems.  The  goal  of  system  identification  is  to  determine  the  parameters  of  an  unknown 
plant  using  measurements  of  the  inputs  to,  and  outputs  from,  the  plant.  A  well-known,  but 
often  conveniently  ignored,  physical  reality  is  that  measurement  entails  uncertainty,  and  the 
measurements  of  the  outputs,  and  possibly  the  inputs,  will  be  corrupted  by  an  unmeasure- 
able  disturbance,  viz.,  measurement  noise.  One  should  also  not  lose  sight  of  the  purpose 
of  system  identification,  which  can  either  be  used  for  merely  determining  the  parameter 
values  of  interest,  or  more  often,  assuming  an  accurate  plant  model,  for  the  purpose  of  de¬ 
signing  a  control  law  for  controlling  the  physical  plant.  Furthermore,  there  are  two  realms 
of  system  identification.  The  first  case  is  a  one  shot,  initial  identification  of  an  unchanging, 
but  unknown,  plant  {e.g.,  for  control  system  design).  The  other  situation  is  the  ongoing 
system  identification  of  a  possibly  time  varying  plant,  in  which  case  one  refers  to  adaptive 
or  reconfigurable  control,  depending  on  whether  the  plant’s  parameters  vary  slowly  or  are 
subject  to  possibly  abrupt  change.  The  latter  is  of  most  interest  in  controls. 

In  principle,  there  are  two  points  of  view  concerning  uncertainty  and  the  control 
of  an  unknown,  or  time  varying,  plant.  One  is  to  rely  on  the  benefits  of  feedback  and 
design  a  robust  feedback  control  system  that  is  valid  for  all  possible  configurations,  and 
the  other  is  to  design  an  adaptive  control  system  that  changes  its  control  based  on  the 
continuously  identified  configuration.  While  (deterministic)  robust  control  can  be  used  to 
guarantee  a  level  of  performance  and  disturbance  rejection,  it  often  leads  to  high  gain  con¬ 
trollers  which  limit  performance  by  saturating  the  actuators.  Moreover,  the  achieved  level 
of  performance  for  a  specific  plant  realization  is  oftentimes  mediocre.  This  can  be  alleviated 
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through  adaptive  control,  where  less  gain  is  required,  viz.,  hard  actuator  saturations  wiU 
be  better  accommodated.  In  addition,  superior  control  performance  will  be  achieved.  The 
need  for  robust  control  is  not  totally  eliminated,  because  the  process  of  system  identification 
is  far  from  perfect.  However,  online  system  identification  reduces  the  level  of  uncertainty, 
so  performance  gains  can  be  realized  through  the  use  of  both  robust  and  adaptive  control. 
Indeed,  the  adaptive/reconfigurable  control  approach  is  most  tempting  in  feedback  control, 
where  no  new  instrumentation  and/or  actuators  beyond  those  already  available  for  robust 
control  are  required.  Obviously,  the  adaptive/reconfigurable  control  compensator  will  be 
nonlinear,  as  opposed  to  the  linear  compensators/controllers  that  are  synthesized  in  the 
current  linear  robust  control  paradigms.  Bearing  in  mind  that  modern  compensators  are, 
in  fact,  algorithms,  and  in  view  of  the  ever  increasing  affordable  computer  power,  the  com¬ 
plexity  of  nonlinear  control  systems  is  no  longer  a  restraint,  and  the  adaptive/reconfigurable 
paradigm  is  increasingly  becoming  a  viable  proposition. 

1.1  A  Survey  of  the  Literature 

System  identification  is  a  broad,  and  hence  a  very  diverse  field,  for  it  lies  at  the  heart 
of  the  empirical  approach  to  mathematical  modeling.  For  example,  a  recent  issue  of  Auto- 
matica  was  devoted  entirely  to  system  identification.  However,  as  the  editor  states  [9],  there 
was  a  distinct  lack  of  papers  on  the  “application  of  parameter  and  system  identification”. 
Indeed,  only  two  of  the  seventeen  papers  even  remotely  qualify.  This  strongly  suggests  that 
system  identification  is  still  in  the  process  of  development.  Two  major  themes  present  in  the 
current  literature  are  the  identification  of  time  varying  parameters  and  the  identification  of 
continuous-time  models  using  discrete- time  observations.  References  [21]  and  [34],  respec¬ 
tively,  are  survey  papers  that  address  the  current  status  of  these  two  important  lines  of 
research.  Additionally,  there  are  many  recent  papers  that  propose  new  methods  of  dealing 
with  the  current  difficulties. 
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In  Ref.  [5]  the  authors  attempt  to  identify  the  physical  parameters  of  a  building,  mod¬ 
eled  by  a  second-order  thermal  network,  for  the  purpose  of  temperature  control.  Derivatives 
and  integrals  of  measured  variables  are  used  to  calculate  the  recursive  least  squares  estimate 
of  combinations  of  the  physical  variables.  There  is  no  real  validation  performed,  because  the 
actual  parameters  are  unknown.  The  data  is  taken  from  a  real  building,  and  the  calculated 
estimates  are  believed  to  be  accurate  because  they  are  “close”  to  those  estimated  in  previous 
work  using  the  same  building.  As  will  be  pointed  out  in  the  sequel,  establishing  a  “correct” 
validation  paradigm  in  system  identification  is  a  crucial  element  of  system  identification 
and  is  a  major  theme  of  this  research. 

In  Ref.  [27],  continuous-time  parameters  are  once  again  estimated,  but  here  the 
observables  are  passed  through  a  modulating  function  before  being  given  to  the  identification 
algorithm,  viz.,  the  input-output  pair  is  filtered.  Careful  input-output  filtering  will  reduce 
the  effects  of  measurement  noise,  and  it  can  result  in  better  estimates.  Unfortunately  the 
method  requires  human  intervention  for  the  design  and  tuning  of  the  filter.  Here,  the  user  is 
required  to  choose  the  “modulating  frequency  index.”  The  estimation  algorithm  is  correctly 
validated  by  running  a  simulation  with  known  parameters  and  measurement  noise,  and  the 
estimated  parameters  are  compared  to  the  known  parameters. 

Another  proposed  method  of  identifying  continuous-time  parameters  is  to  relate  the 
fractional  decomposition  of  the  continuous-time  system  to  several  ^^-transform  polynomials 
[29].  Weighted  least  squares  is  used  to  identify  the  ^-transform  parameters,  and  these 
are  converted  back  to  the  continuous- time  parameter  and  initial  condition  estimates.  The 
continuous-time  model  identification  method  uses  filtering  on  the  input-output  data  pair 
to  estimate  the  models.  A  continuous-time  formulation  of  an  ARMAX  (dynamic)  model  is 
proposed  and  either  recursive  least  squares  or  maximum  likelihood  is  used  for  estimation. 
The  model  is  validated  by  comparing  the  Bode  plots  and  step  and  impulse  responses  of  the 
actual  and  estimated  system. 
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Turning  to  the  time-varying  problem,  in  Ref.  [23]  weighted  least  squares  is  used  to 
examine  the  nonasymptotic  properties  of  finite  memory  identification  as  it  relates  to  time- 
varying  parameters.  This  work  is  based  on  an  ARX  (static)  model  (which  is  not  adequate 
for  dynamic  system  identification),  and  it  discusses  both  forgetting  factors  and  moving 
window  estimation. 

An  improvement  to  the  ad  hoc  forgetting  factor  approach  is  proposed  in  Ref.  [31]. 
In  this  work,  an  online  estimation  of  the  equation  error  covariance  matrix  is  performed 
and  this  estimate  is  used  in  the  recursive  estimation  of  discrete-time  parameters  subject  to 
jumps. 

The  authors  of  Ref.  [35]  are  aware  of  the  problems  associated  with  weighted  recursive 
least  squares  for  time-varying  parameter  estimation,  and  they  propose  a  new  algorithm 
based  on  a  different  parameter  model.  Rather  than  model  the  parameter  as  a  constant, 
they  model  it  as  locally  changing,  and  estimate  the  characteristics  of  the  model.  The 
proposed  algorithm  contains  a  pair  of  user-chosen  tuning  parameters  that  allow  it  to  out¬ 
perform  the  variable  forgetting  factor  recursive  least  squares  algorithm  when  applied  to  a 
discrete-time  ARX  model. 

Model  order  determination  is  also  a  subject  of  interest,  and  the  authors  of  Ref.  [24] 
suggest  an  innovative  method  of  estimating  the  system  order  in  conjunction  with  the  pa¬ 
rameters.  They  augment  and  rearrange  the  parameter  vector  and  regression  matrix  before 
performing  an  LDL^  factorization  of  the  augmented  regression  matrix,  in  which  the  L 
matrix  contains  the  parameter  estimates  for  all  chosen  model  orders,  and  D  gives  the  in¬ 
formation  on  which  model  order  is  most  correct.  Once  again,  the  assumptions  made  are 
consistent  with  an  ARX  model. 

In  Ref.  [18],  the  author  discusses  various  methods  of  parameter  estimation,  including 
maximum  likelihood  estimation,  extended  Kalman  filters,  and  linear  regression  as  they  apply 
to  the  identification  of  correct  models  for  high  performance  aircraft.  Attention  is  given  to 
determination  of  proper  inputs,  correct  model  order,  and  model  validation. 
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Measurement  of  state  derivatives  axe  used  in  Ref.  [6]  to  convert  portions  of  the 
parameter  identification  into  a  “static”  process.  However  low  pass  filters,  with  user  chosen 
roll-off  frequencies,  are  used  elsewhere  in  discretization  equations  for  additional  estimation. 
These  estimates  are  then  used  in  an  adaptive  control  law  for  spacecraft  tracking. 

1.2  Shortcomings  of  the  Current  Identification  Paradigm 

There  are  at  least  three  primary  methods  of  identification.  Maximum  Likelihood 
Estimation  (MLE),  which  is  not  considered  in  this  research,  is  a  sound  method  based  on 
optimization.  It  does,  however,  suffer  from  the  process  of  transcribing  a  difficult  prob¬ 
lem  into  an  “equivalent”  optimization  problem,  which  requires  determination  of  the  global 
minimum.  System  identification  is  one  area  where  the  global  minimum  is  a  requirement, 
because  any  local  minima,  even  those  close  in  value  to  the  global  minimum,  simply  produce 
incorrect  parameter  estimates.  Unfortunately,  global  optimization  is  a  difficult  problem  in 
applied  mathematics. 

Another  classical  route  to  linear  system  parameter  identification  entails  the  estima¬ 
tion  of  the  state  of  an  augmented  and  nonlinear  dynamical  system,  turning  system  identifi¬ 
cation  into  a  nonlinear  filtering  problem.  Hence,  it  would  appear  that  system  identification 
is  in  the  realm  of  the  Extended  Kalman  Filter  (EKE) [25].  In  Extended  Kalman  Filtering, 
a  linearization  is  employed.  However,  when  the  state  estimation  error  becomes  large,  this 
linearization-based  approach  loses  validity  and  the  estimation  algorithm  fails.  One  then 
refers  euphemistically  to  EKF  “divergence.”  In  addition,  the  emphasis  in  Kalman  Filter¬ 
ing  is  on  recursive  algorithms,  and  very  often,  the  complete  measurement  time  history  is 
used.  While  the  recursive  approach  to  estimation  is  most  compatible  with  using  an  ever 
expanding  data  set,  the  latter  has  the  deleterious  effect  of  precluding  the  estimation  of 
time- varying  parameters,  and  in  particular,  the  parameters  subject  to  jumps,  as  is  the  case 
in  systems  subject  to  possible  failures.  Moreover,  EKF’s  must  be  initialized.  When  no  prior 
information  about  the  system’s  states  and  parameters  is  available,  and  the  filter  is  initial- 
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ized  accordingly,  it  might  take  a  long  time  for  the  erroneous  information  to  be  “washed 
out.”  Furthermore,  proper  identification  requires  a  data  window  of  minimal  physical  (tem¬ 
poral)  length.  Estimates  produced  using  windows  shorter  than  this  minimal  length,  even 
though  high  sampling  rates  can  produce  many  data  points,  are  useless.  In  conclusion,  this 
identification  method  suffers  from  the  well  known  deficiencies  of  EKF’s. 

Therefore,  in  the  control  community,  linear  regression-based  approaches  to  the  iden¬ 
tification  of  linear  system  parameters  are  sought.  These  are  also  referred  to  as  “linear 
prediction”  algorithms.  Here,  the  linear  structure  of  the  dynamics  is  directly  exploited  and 
the  system  parameters  only  (without  the  states)  are  estimated.  This  main  line  of  research 
in  system  identification  is  based  on  the  statistical  method  of  linear  regression.  This  prag¬ 
matic  approach  to  system  identification  is  also  the  subject  of  this  research.  The  interest  in 
regression-based  system  identification  arises  from  the  relatively  easier  (compared  to  MLE) 
computational  burden  and  its  relatively  reduced  sensitivity  to,  or  possibly  no  need  for,  an 
initial  parameter  guess.  The  latter  holds  the  promise  of  autonomous  operation,  without  the 
need  for  human  intervention,  also  known  as  “tuning”. 

The  main  theory  behind  regression-based  identification  is  least  squares.  There  are 
many  variations  on  this  theme:  weighted  least  squares,  generalized  least  squares,  instru¬ 
mental  variables,  minimum  variance,  etc.  Each  suffers  from  its  own  difficulties  when  the 
estimation  problem  is  not  properly  formulated.  Unfortunately,  there  is  a  glaring  weakness 
in  the  whole  of  system  identification  today:  the  lack  of  proper  treatment  of  measurement 
noise.  Least  Squares  (LS),  which  is  widely  used  in  statistics,  is  a  sound  method  as  long  as 
the  assumptions  upon  which  it  is  based  are  met.  Unfortunately,  as  will  be  discussed  in  the 
sequel,  the  identification  of  dynamic  systems  with  measurement  noise  negates  an  important 
noncorrelation  assumption  in  LS. 

The  problem  is  further  aggravated  by  the  lack  of  proper  treatment  of  dynamical 
system  identification  with  measurement  noise  in  some  of  the  most  cited  textbooks.  When 
regression  based  identification  is  discussed  in  Refs.  [12],  [14],  and  [20],  measurement  noise 
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is  barely  mentioned  as  the  discussion  moves  from  least  squares  theory  to  dynamical  system 
identification.  Unfortunately,  the  presence  of  measurement  noise  in  dynamical  systems 
produces  a  correlation  that  biases  the  LS  estimate.  The  LS  variations  mentioned  above  were 
formulated  to  address  the  correlation  problem,  and  it  could  be  assumed  that  measurement 
noise  is  just  lumped  together  with  other  correlated  noises.  However,  this  is  not  a  judicious 
course  of  action,  because  measurement  noise  introduces  a  specific  correlation  that  can  be 
modeled  and  successfully  addressed. 

Another  misplaced  emphasis  in  the  literature  is  on  recursive  estimation.  The  benefits 
of  this  methodology  are  somewhat  exaggerated.  It  is  true  that  the  ongoing  estimation  of 
a  parameter  vector  could  possibly  be  accomplished  faster,  and  with  decreased  memory 
requirements,  with  a  recursive  method,  but  care  must  be  taken  in  the  application  of  such 
methods.  In  recursive  estimation,  an  initial  guess  is  required  as  to  the  value  of  the  parameter 
and  the  confidence  in  the  guess.  If  there  is  no  a  priori  knowledge  of  the  parameter,  typically 
a  random  guess  is  made,  and  the  confidence  in  it  is  set  very  low.  This  may  seem  like 
a  proper  thing  to  do,  but  poor  initial  guesses  can  take  a  long  time  to  wash  out  of  the 
recursively  calculated  estimate.  Another  fallacy  of  the  recursive  methods  is  that  one  obtains 
an  estimate  after  the  first  measurement,  but  this  estimate  is,  in  fact,  meaningless.  A 
certain  number  of  measurements,  determined  by  a  minimal  p/ij/s?ca/ identification  interval,  is 
required  to  obtain  a  valid  estimate,  and  that  number  does  not  decrease  when  using  recursive 
methods.  In  fact,  it  could  increase  because  of  the  poor  initial  guess.  At  the  same  time,  the 
recursive  formulation  of  an  identification  algorithm  oftentimes  obscures  important  aspects 
of  its  operation.  This  aspect  of  useful  estimates  will  also  be  addressed  in  the  proposed 
research. 

The  identification  of  time- varying  parameters  is  often  based  on  recursive  estimation. 
The  fact  that  the  parameters  vary  wiU  produce  an  incorrect  estimate  as  time  progresses. 
This  is  because  the  method  of  recursion,  by  its  nature,  remembers  all  the  information.  An 
ad  hoc  method  of  “losing”  the  initial  information  is  to  introduce  a  forgetting  factor  into 
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the  recursion.  If  it  is  true  that  the  parameter  is  varying,  a  better  method  of  estimation 
would  be  a  moving  identification  window.  The  latter  method  is  more  amenable  to  analysis, 
and  the  window  length  can  be  more  directly  determined  based  on  information  about  the 
physical  parameter’s  drift. 

There  are  two  methods  of  continuous- time  system  parameter  identification  [34].  One 
is  to  estimate  the  parameters  of  a  discretized  (in  time)  system  using  weU-known  regression 
techniques,  and  then  convert  the  discrete-time  system  back  into  a  continuous-time  system 
using  a  (bilinear)  transformation.  A  problem  arises  here,  which  is  addressed  in  the  disserta¬ 
tion,  concerning  the  behavior  of  the  sampled  input  signal.  After  all,  the  continuous  physical 
plant  is  subject  to  a  continuous  input  signal,  whereas  only  its  sampled  values  are  available 
to  the  identification  algorithm.  The  other  method  of  continuous-time  system  identification 
entails  modeling  the  continuous  parameters  directly  and  then  generating  certain  “measures” 
from  the  discretely  observed  input  and  output  that  can  be  related  to  the  parameters  in  a 
regression-based  manner.  A  troublesome  aspect  of  this  method  is  the  generation  of  these 
“measures.”  It  requires  human  intervention  to  select  certain  parameters  involved.  Gener¬ 
ating  measures  is  equivalent  to  “properly”  filtering  the  input /output  pair.  A  truly  rigorous 
treatment  of  the  continuous-time  case,  with  continuous-time  measurements,  resides  in  the 
realm  of  the  difficult  mathematical  theory  of  nonlinear  filtering. 

1.3  Research  Objectives 

One  of  the  main  themes  of  this  dissertation  is  that  system  identification  is  not  merely 
an  algorithm  which  can  be  universally  and  blindly  applied  with  no  prior  knowledge.  Thus, 
the  main  tenant  of  this  work  is  that  system  identification  is,  in  fact,  a  process  where  the 
following  factors  play  a  major  role: 

1.  The  recognition  that  measurement  noise  is  ubiquitous  in  physical  systems 

2.  The  availability  of  prior  information  about  the  possible  system  order  and  parameter 
values  and  confidence  bounds  on  the  latter 
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3.  Relevant  dynamics  and  the  bandwidth  of  interest 


4.  The  presence  of  parasitic  dynamics,  both  high  and  low  frequency 

5.  The  way  in  which  these  parasitic  dynamics  interact  with  the  dynamics  of  interest 

6.  Excitation  issues,  including: 

•  Physical  duration  of  the  identification  interval  T 

•  Sampling  rate  1/AT  or  number  of  measurements  m  =  T/AT 

•  The  number  of  identified  parameters  must  be  small  as  compared  to  m  (use  lower 
order  models  if  possible) 

•  Signal  to  Noise  Ratio  (SNR) 

•  Frequency  content  of  the  input  signal 

The  main  thrust  of  this  research  is  to  develop  a  “correct”  and  comprehensive  system 
identification  paradigm  which  addresses  the  items  above.  To  this  end,  proper  stochastic 
modeling  of  the  measurement  noise  is  of  paramount  importance.  Initially,  the  theory  of 
minimum  variance  estimation  is  applied  in  an  innovative  way,  based  on  the  correct  noise 
model,  and  examined  in  the  face  of  input  and  output  measurement  noise.  Further  subjects 
of  interest  are  the  application  of  this  concept  to: 

1.  Phasor-based  system  identification  of  continuous-time  systems 

2.  Static  identification  of  continuous-time  systems 

3.  Time  domain  identification  of  discrete-time  systems 

Although  the  estimation  is  performed  with  discrete,  noise-corrupted  measurements,  the 
primary  goal  of  this  research  is  to  identify  the  physical  parameters  of  a  plant.  For  that,  the 
continuous-time  model  is  required. 

Moreover,  in  this  work,  the  emphasis  is  on  system  identification  under  the  following 
restrictions: 
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•  Small  sample;  Proof  of  asymptotic  results  (as  m  oo)  is  quite  useless.  For  one  thing, 
the  noise  is  then  washed  out.  Furthermore,  in  practice,  the  small  sample  identification 
horizon  is  typically  too  short  for  the  asymptotic  results  to  have  any  bearing  on  the 
problem. 

•  Real  time  operation;  The  purpose  of  the  proposed  identification  is  for  adaptive  control. 
Therefore,  regression-based  system  identification,  which  can  be  adapted  to  online 
operation,  is  used,  rather  than  maximum  likelihood  estimation. 

•  Autonomous  operation;  Identification  will,  ideally,  be  performed  without  human  in¬ 
tervention  and  tuning. 

•  Measurement  noise;  It  is  understood  that  this  is  a  reality  of  physical  system  identifi¬ 
cation. 

The  primary  consequence  of  the  restrictions  above  is  that  the  number  of  identified  parame¬ 
ters  must  be  kept  small.  It  is  impossible  to  identify  large  numbers  of  parameters  accurately 
with  a  small  sample  and  significant  measurement  noise. 

Finally,  above  and  beyond  the  system  identification  algorithm  development  work, 
attention  is  given  to  the  validation  process,  especially  in  the  realistic  case  where  one  operates 
on  real  data. 

1.4  Methodology 

To  achieve  these  objectives,  a  novel  system  identification  algorithm  is  developed  and 
used  in  this  work.  Here,  the  algorithm  is  coined  Generalized  Minimum  Variance  (GMV), 
because  of  similarities  to  the  Generalized  Least  Squares  algorithm,  and  the  use  of  Mini¬ 
mum  Variance  estimation  in  the  solution.  Introduced  in  this  work  is  the  concept  that  the 
parameter  estimate  provided  by  the  GMV  algorithm  is  a  fixed  point  of  a  nonlinear  map¬ 
ping  derived  from  Minimum  Variance  estimation.  The  existence  of  a  fixed  point  is  proven. 
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and  the  convergence  properties  of  the  algorithm  are  examined  for  difficult  identification 
problems  with  low  signal-to- noise  ratios. 

The  novel  GMV  algorithm  was  first  used  in  Ref.  [4]  for  determining  the  optimal 
inputs  for  the  identification  of  discrete-time  plants  with  output  measurement  noise  only. 
This  noise  scenario  is  examined  in  this  work,  and  expanded  to  include  input  measurement 
noise  as  well. 

In  addition  to  the  discrete-time  identification  problem,  the  GMV  algorithm  is  applied 
to  the  frequency  domain  identification  of  a  continuous-time  plant,  using  three  different 
noise  scenarios.  Also  examined  in  this  work  is  the  identification  of  the  parameters  of  a 
physical  system  using  samples  of  the  continuous-time  inputs  and  output.  In  each  case,  the 
performance  of  the  GMV  algorithm  is  examined  under  decreasing  signal-to-noise  ratios. 
Careful  and  weU-documented  tests  are  performed,  and  the  results  are  compared  to  several, 
more  common,  system  identification  algorithms. 

1.5  Organization  of  the  Dissertation 

Chapter  II  focuses  on  the  idea  of  linear  regression  (or  linear  prediction)  for  system 
identification.  The  basic  theory  of  Least  Squares  is  discussed,  along  with  the  introduction 
of  correlated  equation  error.  The  latter  is  a  result  of  measurement  error  incurred  during 
work  with  dynamical  systems,  and  is  revealed  by  proper  modeling  of  the  measurement 
process.  Proper  treatment  of  the  equation  error  correlation  is  then  developed,  and  the 
theory  behind  Minimum  Variance  identification  is  applied  to  the  problem.  This  leads  to 
the  GMV  algorithm. 

Chapter  III  contains  the  application  of  the  GMV  algorithm  to  the  identification  of 
continuous-time  transfer  functions  using  frequency  domain  measurements.  Several  different 
measurement  models  are  considered.  These  include  a  basic,  but  unlikely,  model  used  for 
comparison  purposes,  and  a  realistic  measurement  model  that  reflects  the  capabilities  of 
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current  frequency  analyzers.  The  GMV  identification  results  are  compared  to  the  results  of 
other  more  common  identification  methods. 

Chapter  IV  is  a  short  discussion  of  the  effects  of  unmodeled  dynamics  on  the  fre¬ 
quency  domain-based  identification  method.  For  example,  current  flight  control  system 
design  is  often  performed  on  low  order  models,  when  in  fact  the  true  plant  contains  both 
unmodeled  low  and  high  order  dynamics. 

Chapter  V  contains  the  application  of  the  GMV  algorithm  to  the  identification  of 
discrete-time  dynamical  systems  using  discrete-time  measurements.  The  assumption  here  is 
that  the  plant  is  being  controlled  with  a  zero-order  hold,  so  the  continuous-time  parameters 
can  be  obtained  simply  by  performing  an  inverse  zero-order  hold  transformation.  The 
GMV  algorithm  results  are  compared  at  each  step  to  those  of  alternative,  more  common, 
identification  schemes.  The  effects  of  initial  transients,  input  and  output  noise,  and  sampling 
rate  are  aU  investigated  here. 

Chapter  VI  is  a  discussion  of  the  identification  of  a  continuous-time  system’s  physical 
parameters  using  sampled  measurements  of  the  system’s  input  and  output.  This  work  is 
based  on  a  competition  that  originated  in  Italy.  The  basic  dynamical  system  being  identified 
is  based  on  the  thermal  characteristics  of  buildings,  but  it  is  basically  modeled  as  a  simple 
second-order  electrical  circuit. 

Finally,  chapter  VII  is  a  summary  and  discussion  of  the  results  of  the  research.  This 
is  a  rich  field  of  work,  and  this  research  should  only  be  a  beginning. 
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II.  Linear  Regression 


In  the  system  identification  paradigm,  there  is  an  observable  Z,  that  is  related  to  an 
unknown,  but  desired,  parameter  vector  0  . .  .OnY •  If  ^  is  constant,  and  that  relation 

is  linear,  then  a  number  of  observations  (m)  can  be  taken 

Zk  =  Hkl^l  +  Hk2^2 - \- Hkn^n  *=l,2,...,m  (2-1) 

where  (2.1)  is  referred  to  as  a  regression  function  [14].  These  observations  can  then  be 
arranged  in  a  matrix  form  linear  regression 
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Here  H  is  referred  to  as  the  regression  matrix.  There  are  many  methods  to  estimate  0  based 
on  the  observations,  and  most  are  related  to  Least  Squares  (LS). 

2. 1  Least  Squares  and  Weighted  Least  Squares 

The  estimation  method  of  least  squares  was  originated  by  Gauss  in  the  early  19th 
century,  when  he  used  it  to  predict  the  orbits  of  planets.  Since  then,  it  has  become  a 
commonly  used  tool  in  parameter  estimation.  Its  popularity  among  engineers  stems  from  the 
fact  that  it  is  easier  to  understand  than  some  other  methods,  such  as  maximum  likelihood, 
and  that  it  does  not  require  a  knowledge  of  mathematical  statistics  [14].  “Deterministic” 
optimization  approaches  abound.  Also,  for  certain  problems  that  are  formulated  properly, 
the  least  squares  method  yields  estimates  which  are  consistent,  unbiased,  and  efficient. 

If  there  was  no  uncertainty  in  the  observations  {Zk  and  Hki  obtainable),  only  m  =  n 
linearly  independent  observations  would  be  needed  to  calculate  the  parameters  uniquely. 
However,  “all  our  measurements  and  observations  are  nothing  more  than  approximations 
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to  the  truth  [10],”  so  there  is  uncertainty  present,  and  this  requires  one  to  make  m  >  n 
observations.  At  the  very  least,  there  is  uncertainty  about  the  Zfe’s,  he.,  one  can  obtain  an 
uncertain  Zk  =  Zk  +  v^.  Therefore, 

=  HklOi  +  irjfc2^2  + - h  Hkn^n  + 

which  implies  that 

z  =  H0  +  e  where  e  =  [ui . . .  Vm]'^  and  z  =  [zi . . .  Zm]'^ 

In  least  squares  theory,  the  objective  is  to  find  the  estimate  0  which  minimizes  the  square 
of  the  equation  error  vector  e{0)  =  z  —  HO,  he.,  [14] 

=  argmine(0)^e(0)  =  H^z  (2-2) 

A  basic  assumption  in  the  above  formulation  is  that  the  uncertainty  in  each  of  the 
observations  is  the  same.  If  it  is  not,  then  a  weighting  matrix  W  can  be  introduced  to 
weight  each  observation  differently,  and  the  Weighted  Least  Squares  (WLS)  estimate  is 
then  given  by 

Owls  =  arg  min  e{0fWe{e)  =  (h^Wh)  H^Wz  (2.3) 

If  some  prior  knowledge  of  e  exists,  specifically  that  e  is  a  stationary  random  vector  with 
£[e]  =  0  and  E[€e^]  =  R,  where  E[-]  denotes  statistical  expectation,  then 

0MV  =  (h^R"^h)”^  H^R-^z  (2.4) 

is  the  minimum  Best  Linear  Unbiased  Estimate  (BLUE)  [12],  or  the  estimate  which  has 
the  minimum  variance  out  of  all  linear  unbiased  estimates.  Therefore,  it  is  referred  to  as 
the  Minimum  Variance  (MV)  estimate.  For  the  special  case  when  the  Vk,  k  =  1, . .  .,m, 
are  identically  distributed  and  independent  (i.i.d.)  with  variance  cr^,  the  covariance  of  the 
estimation  error  is  a  scaled  identity  matrix  (R  =  cr^I),  and  the  estimate  0mv  =  0LS- 
Alternatively,  if  one  assumes  that  Vi  are  i.i.d.  and  Gaussian,  then  0mv  is  the  maximum 
likelihood  estimate  of  0  [12]. 
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The  primary  drawback  of  LS  estimation  is  the  requirement  that  R  =  cr^I.  If  it  is  not, 
as  in  the  case  when  measurement  noise  is  present  and  dynamical  systems  are  considered, 
i.e.,  the  traditional  control  situation,  then  6]^s  will  be  biased.  The  WLS  method  eliminates 
the  bias,  but  R  is  often  not  known  in  practice.  One  of  the  main  thrusts  in  this  research  is 
the  correct  modeling  of  R,  yielding  parameter  estimates  with  smaller  bias. 

2.2  Correlation  Caused  by  Dynamical  Systems 

In  this  work,  the  primary  interest  lies  in  the  identification  of  dynamical  systems,  and 
in  this  type  of  problem,  not  only  are  the  Zfc’s  unobtainable,  but  also  the  Rfci’s.  Now,  the 
obtainable  quantity  is  hki  —  Hki  +  Wkii  and  the  pertinent  equation  is 

Zk  =  hklO-i  +  +  •  •  •  +  hknOn  +  - WkJn  (2.5) 

Here,  6  =  [ei . .  .6^]^  where  Ck  =■  Vk  —  ^i'^ki-  It  is  known  that,  as  m  — >■  oo,  the  least 

squares  estimate  converges  as  [12] 

Ols  —  O^-  (2.6) 

If  H  and  e  are  uncorrelated,  then  ^[H^e]  =  H^i?[€]  =  H^O  =  0,  but  for  the  system  given 
in  Eq.  (2.5),  the  fcth  element  of  E\R'^e\  is  given  by 

hik 

Li=i 

Even  if  E[wkiWkj]  =  substituting  for  hik  above  yields 

j  =  ma^9k  0 

and  this  correlation  causes  the  LS-based  estimate  to  be  biased.  Hence,  when  measurement 
noise-corrupted  dynamical  models  are  used,  it  is  important  to  recognize  the  correlation 
inherent  in  the  linear  regression’s  equation  error.  It  is  therefore  wise  to  estimate  using  a 
method  which  incorporates  and  models  the  equation  error  covariance  information. 


m 
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2.8  Dealing  with  Correlation 


There  are  several  methods  in  the  literature  today  that  were  developed  to  overcome 
the  correlation  introduced  by  the  dynamical  model,  thereby  reducing  the  bias  in  the  LS 
estimate.  The  two  mentioned  here  are  the  Generalized  Least  Squares  (GLS)  and  the  In¬ 
strumental  Variable  (IV)  methods.  However,  these  methods  do  not  specifically  address  the 
particular  correlation  introduced  by  measurement  noise  in  dynamical  systems. 

2.3.1  Generalized  Minimum  Variance.  To  address  the  issue  of  measurement  noise 
properly,  the  following  procedure  is  offered.  The  available  information  can  be  arranged  in 
a  linear  regression  equation  given  by 

z  =  H0 -I- €  €  =  T(0)v,  v  =  A/^(0,R,;),  R„  >  0  (2.7) 

To  use  insight  gained  from  the  minimum  variance  identification  method,  it  is  necessary  to 
obtain  the  equation  error  covariance  R,  where 

R  =  E[ee'^]  =  TR,,T^  (2.8) 

Unfortunately,  R  is  not  known  a  priori,  because  in  addition  to  the  dependence  on  the 
given  sensor’s  measurement  error  cr„,  it  is  a  function  of  the  (as  yet  unknown)  coefficients 
of  the  system’s  transfer  function,  i.e.,  R  =  R(0).  Two  related,  but  different,  estimates 
can  be  derived  from  this  minimum  variance  based  approach.  First,  one  could  minimize  the 
associated  cost  function  to  obtain  the  estimate: 

6  —  argmine(0)^€(0)  =  (z  —  H0)^R“^(0)  (z  -  H0)  (2.9) 

0 

Because  of  the  dependence  of  R  on  6,  this  leads  to  a  complicated  numerical  search  for 
a  global  minimum.  In  global  nonlinear  searches  such  as  this,  the  appearance  of  local  minima 
is  virtually  guaranteed  as  the  noise  levels  increase.  There  is  a  problem  with  local  minima 
in  system  identification  problems  that  tends  not  to  be  a  problem  in  other  minimization 
problems.  While  a  local  minimum  in  a  general  minimization  problem  probably  produces  a 
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solution  with  a  lower  cost,  which  may  be  acceptable  as  a  solution,  the  local  minimum  in 
a  system  identification  problem  can  be  very  far  from  the  required  global  minimum  and  is 
probably  incorrect.  Additionally,  in  general  problems,  a  local  minimum  with  a  cost  that  is 
very  close  to  the  global  minimum  cost  probably  produces  an  acceptable  solution.  However, 
the  fact  that  the  cost  at  a  local  minimum  in  a  system  identification  problem  is  close  to 
the  global  minimum  cost  has  no  bearing  on  the  “closeness”  of  the  estimate  to  the  correct 
solution.  It  most  likely  results  in  a  completely  different,  and  incorrect,  estimate.  Therefore, 
only  the  global  minimum  is  an  acceptable  solution  in  system  identification. 

In  an  attempt  to  avoid  this  complication,  Eq.  (2.4)  is  used  to  obtain  the  second 
possible  derivation  of  the  estimate.  The  Generalized  Minimum  Variance  (GMV)  estimate 
is  given  by  the  point  Ogmv  such  that 

Bgmv  =  (h^R-H«gw)h)~'  H^R-1(0gmv)z  (2.10) 

There  are  many  different  ways  of  searching  for  fixed  points  like  Ogmv,  but  the  following 
simple  iterative  algorithm  has  been  effective  in  finding  the  correct  fixed  point  for  signal-to- 
noise  ratios  approaching  0  dB.  This  algorithm  is  motivated  by  the  iteration  for  fixed  points 
of  contraction  mappings,  for  which  the  existence  of  a  fixed  point  is  guaranteed  [1]. 

Step  1  -  Set  i  =  0  and  calculate  an  initial  parameter  estimate  using  LS. 

=  00  =  H^z 

Step  2  -  Calculate  R(08)  using  Eq.  (2.8). 

Step  3  -  Calculate  9i+i  via  Eq.  (2.4). 

di+i  =  (H^R-i(0i)H)'^H^R-H0i)z 

Step  4  -  If  II  Oi^i  —  0j  II  is  less  than  some  acceptable  value,  proceed  to  step  5.  Otherwise, 
increment  i  and  return  to  step  2. 

Step  5  -  Set  Ogmv  —  ^j+i- 
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step  6  -  The  error  covariance  of  the  estimate  Ogmv  is  then  given  by 

(h^R-1(0gmv)h)"'  (2.11) 

A  fixed  point  for  Eq.  (2.10)  does  exist,  at  least  for  the  problems  of  linear  dynamical 
systems,  as  is  shown  in  Section  3.3.  Also,  the  algorithm  has  converged  within  the.  numerical 
limits  of  Matlab  for  aU  problems  examined  thus  far.  The  number  of  iterations  required  is 
quite  small  for  small  noise  levels,  but  does  increase  as  the  noise  level  increases. 

Problems  do  arise  as  the  noise  levels  increase,  as  is  common  in  system  identification 
problems  at  low  signal-to-noise  ratios.  In  this  work,  it  is  common  for  one  other  fixed  point 
to  appear  at  higher  noise  levels,  and  possible  methods  of  dealing  with  this  are  discussed  in 
Sections  3.6.3  and  5.4.  However,  the  other  methods  of  identification  used  for  comparison  in 
this  work  tend  to  suffer  at  these  noise  levels  as  weU. 

2.3.2  Other  Methods  of  Correlation  Compensation.  There  are  many  methods  of 
compensating  for  correlated  measurement  noise,  as  discussed  in  [30],  but  the  two  used  for 
comparison  in  this  work  are  the  Instrumental  Variable  (IV)  method  and  the  Generalized 
Least  Squares  (GLS)  method.  These  two  methods,  or  variations  thereof,  are  common  in 
the  literature.  One  problem  with  each  of  them  is  that  the  derivation  tends  to  be  specific  to 
the  given  problem,  so  the  particular  method  must  be  discussed  in  context,  as  is  done  in  the 
sequel. 

Briefly  put,  however,  the  IV  method  uses  the  fact  that  the  Zfc’s  and  Lffci’s  are  actually 
functions  of  some  fixed,  but  unknown  B.  During  the  IV  identification  method,  the  current 
estimate  is  used  to  produce  a  set  of  estimated  which  are  then  used  in  an  iterative 

scheme  where 

0i+i  =  H(Ii)^z 

Similarly,  the  GLS  method  exploits  the  form  of  e  in  Eq.  (2.7).  The  contents  of  the 
z  and  H  matrices  are  filtered  through  some  form  of  T~^,  and  the  LS  estimate  is  calculated 
using  these  filtered  results.  The  filtering  tends  to  “whiten”  the  problem,  reducing  the 
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correlation  and  causing  the  LS  formulation  to  be  more  correct.  This  is  also  an  iterative 
procedure,  as  T  is  a  function  of  6,  so  the  contents  of  z  and  H  must  be  re-filtered  at  each 
step. 

2.4  Conclusions 

The  correlation  introduced  by  dynamical  systems  and  measurement  noise  effectively 
eliminates  the  LS  estimate  as  a  valid  result.  Some  method  of  dealing  with  this  correlation 
must  be  found  in  order  to  obtain  an  accurate  estimate.  There  are  many  current  methods  of 
dealing  with  the  correlation,  but  none  seem  to  specifically  address  the  particular  correlation 
introduced  by  measurement  noise. 

The  stochastic  modeling  and  GMV  identification  concept  described  in  this  chapter 
are  not  limited  to  time-domain  dynamical  systems.  They  can  also  be  applied  to  any  problem 
where  certain  measurements  are  correlated  with  one  another.  An  example  of  this  is  the 
frequency  domain  identification  of  a  continuous- time  plant,  which  is  discussed  in  the  next 
chapter. 
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III.  Phasor-Based  Identification  of  a  Continuous-Time  Dynamical  System 


In  this  chapter,  a  phasor  approach  to  system  identification  is  discussed,  viz.,  the 
experimental  data  consists  of  a  finite  number  of  sensor  noise  corrupted  point  frequency 
response  measurements  of  the  unknown  plant. 

The  proposed  identification  paradigm  is  in  line  with  currently  available  instrumen¬ 
tation,  e.g.,  frequency  analyzers  [33].  Gaussian  measurement  noise  statistics  are  assumed. 
The  emphasis  here  is  on  a  stochastic  analysis  of  the  identification  problem  in  the  frequency 
domain.  Parameter  estimation  algorithms  are  developed  and  validated  in  simulations  which 
include  measurement  noise.  Careful  modeling  of  the  stochastic  estimation  problem  renders 
an  efficient  system  identification  algorithm. 

Moreover,  the  proposed  system  identification  algorithm  yields  an  identified  model 
and  an  estimate  of  the  model  uncertainty.  The  error  is  expressed  in  terms  of  the  uncertainty 
of  the  coefficients  of  the  plant’s  transfer  function,  the  latter  being  easily  transformable 
into  an  expression  of  uncertainty  about  the  physical  parameters  of  the  plant.  In  other 
words,  the  proposed  system  identification  algorithm  directly  yields  an  identified  model  and 
its  structured  uncertainty,  as  opposed  to  plant  uncertainty  gauged  with  the  Hoo  norm. 
Furthermore,  the  Hoo  estimation  error  can  be  calculated  using  optimization  methods,  while 
Kharitonov  [17]  type  results  can  be  invoked  to  ascertain  robust  stability.  Obviously,  the 
measurement  error  is  somewhat  indicative  of  the  Hoo  estimation  error. 

The  phasor  approach  is  developed  in  Section  3.1.  A  stochastic  analysis  of  the  sys¬ 
tem  identification  problem  with  measurement  noise  is  performed  in  Section  3.2.  Section  3.3 
contains  a  proof  that  at  least  one  fixed  point  exists  for  this  formulation.  Section  3.4  ex¬ 
plains  how  to  calculate  a  needed  covariance  matrix  for  different  assumptions  about  the 
measurement  model.  A  brief  discussion  of  the  other  algorithms  being  compared  is  con¬ 
tained  in  Section  3.5.  Sections  3.6  contains  the  results  of  applying  the  GMV  algorithm  to 
a  second-order  system,  as  well  as  the  comparisons  to  the  other  algorithms.  The  system 
used  is  representative  of  the  short  period  dynamics  of  an  aircraft.  The  performance  of  the 


3-1 


different  system  identification  algorithms  is  briefly  discussed  in  Section  3.7,  and  concluding 
remarks  are  made  in  Section  3.8. 


3. 1  The  Frequency  Domain  Based  Phasor  Approach 

Consider  the  stable  frequency  domain  transfer  function,  where  the  order  of  the  nu¬ 
merator  and  denominator  are  known: 

T(s)  ~  ^  ^  + - 1~  bn-lS  +  bn 

u(s)  — - On-l-S  ~  O-n 

Here,  the  parameters  . . .  a„  and  h\. .  .b^  are  unknown,  the  a\. .  .a^  coefficients  are  neg¬ 
ative  for  a  stable  system,  and  the  b\...  6g_i  coefficients  equal  zero  when  the  order  of  the 
numerator  is  q  less  than  the  order  of  the  denominator.  If  these  parameters  are  equal  to  zero, 
it  should  be  taken  into  account  when  the  problem  is  developed.  If  one  attempts  to  blindly 
identify  all  the  parameters,  assuming  that  the  zero  value  parameters  wiU  just  identify  to 
zero,  the  results  will  suffer. 

By  applying  an  input  to  this  system  of  the  form 

Uk{t)  =  COs{(jJkt) 

and  waiting  for  steady-state  to  be  achieved,  the  steady-state  output  is 
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Through  algebraic  manipulations,  one  obtains  the  linear  system; 


L  J 

(3.3) 

As  can  be  seen,  each  sinusoidal  test  input  produces  two  equations  in  the  2n  un¬ 
knowns.  Therefore,  n  distinct  sinusoids  are  needed  to  produce  the  2n  equations  required 
for  the  determination  of  the  parameter  vector  0  =  [oi  . . .  6i  ...  G  7^^". 

8.2  Measurement  Noise 

In  reality,  no  physical  measurements  can  be  made  without  some  sort  of  measurement 
noise  present,  so  the  true  Ak  and  Bk  are  not  available.  Rather,  there  is  the  noise  corrupted 
Ak^  and  Bk„ ,  where  the  measured  quantities  are 

Akm  =  Afc  +  VA^ ,  va  =  ^"(0,  or  a) 

Bkm  =  Bk  A  UBfc,  vb  =  Af{Q,  cr|) 

Here,  the  random  variables  va,,  and  vb^  are  the  measurement  noises  in  the  fcth  experiment. 


fc  =  1, 2, . . . ,  m.  Therefore  in  reality,  Eq.  (3.3)  is 


“ 

n 

n 

Wk 

{Skm—^Bk) 
- ^ - 

-(^km-^Ak) 

JL 

Wk 

0  4  ••• 

0  = 

-{Bkm 

-^Bk) 

OJh 

••  0 

4  0  ••• 

-^Ak) 

*■ 

-* 

(3.4) 

This  equation  can  now  be  written  as  a  regression  function  with 
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where  the  2x2  Rfc  submatrix  (k  =  1, 2, . . . ,  m)  is  of  the  form 


Rfc  =  Ejefeel’}  =  E^TkVkvlTl'^ 


=  T^RabTl  where  R^b  = 


^AB 
OAB  0's 


Note  that 


R-^  = 


R7 


0  •  •  •  R: 


,  H^R-^H  =  Y.  R-^z  =  Y 

fc=i  fe=i 


so 


-1 


i  nt,  \  -  nt> 

OqMV  -  I  Y  HrRfc^(^GMv)Hfc  j  Y  HfcRfc^(0aMv)ZA: 

\A:=1  /  k=\ 

Pgmv  =  ^KlRk\eGMv)nk^ 


(3.7) 


3.S  Existence  of  Fixed  Point 


The  primary  problem  in  solving  for  the  GMV  estimate  is  that  R  is  a  nonlinear 
function  of  0.  Therefore,  define  /  :  7^^"  — >■  7?.^", 

f{0)  =  [h^R-1(0)h]'^  H^R-1(6>)z  (3.8) 

and  use  an  iteration  to  search  for  a  fixed  point  in  7^^"  where  9  =  f{9).  In  searching  for 
a  fixed  point,  it  would  be  nice  to  know  that  a  fixed  point  actually  exists.  If  it  could  be 
determined  that  f{0)  is  a  continuous,  bounded  function,  then  it  would  be  easy  to  show  that 
at  least  one  fixed  point  exists. 


3.3.1  Singularities  and  Continuity.  The  first  step  in  showing  that  f{6)  is  bounded 
is  to  show  that  it  is  a  continuous  function,  or  that  it  has  only  removable  singularities.  To  do 
this,  first  return  to  the  equation  error  noise  vector  in  Eq.  (3.6).  Because  sums,  products,  and 
inverses  of  rational  functions  are  rational  functions,  and  Tfc  is  a  matrix  of  rational  functions 
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in  the  parameters,  then  f{6)  is  a  rational  function.  Therefore,  f{6)  is  continuous  everywhere 
it  is  defined.  However,  f(0)  is  not  defined  everywhere,  but  has  singularities  where  R  is  not 
invertible.  Since  R^b  is,  by  definition,  positive  definite,  the  singularities  occur  when  the 
determinant  of  any  particular  Tfc  matrix  evaluates  to  zero.  Since  the  transformation  matrix 
Tfc  is  affine  in  the  denominator  parameters,  it  can  be  written  as 


—  Ao  +  ^2  Ajfcttj  (3-9) 

j=i 

From  Eq.  (3.6),  it  can  be  seen  that  T  always  has  the  form 

a  b 

T  = 

—b  a 

This  implies  that  det(T)  =  +  b^,  and  the  only  way  that  this  can  be  zero  is  if  a  and  b 

are  zero,  i.e.,  T  =  0,  For  the  case  n  =  1,  it  is  impossible  for  Tk  to  equal  0.  For  n  =  2, 
the  only  way  T*.  can  equal  0  is  if  oi  =  0  and  02  =  This  leads  to  m  distinct  singular 
points  of  f{6);  one  for  each  measurement  frequency.  For  n  =  3,  if  02  =  and  03  =  w|ai, 
then  Tfc  =  0.  This  implies  that  there  are  m  distinct  lines  of  singularities  in  0  space,  or 
m  distinct  singular  manifolds  of  dimension  1.  For  general  n,  the  dimension  of  the  singular 
manifolds  is  n  —  2. 

To  get  an  idea  of  the  nature  of  these  singularities,  an  n  =  2  case  is  considered.  Three 
measurements  (w  =  1, 2, 10)  are  taken  of  a  second-order  plant,  and  f(0)  is  plotted  along  the 
path  Cl  =  0.  The  results  in  Fig.  3.1  appear  to  show  that  the  singularities  are  removable. 
Indeed,  this  seems  to  be  the  case  for  higher  dimension  problems  as  well.  If  f{0)  is  evaluated 
at  0  -b  e,  where  0  is  a  singular  point  and  e  is  some  small  vector,  the  results  are  close  for  any 
€,  implying  that  the  singularity  is  removable  as  well.  For  the  singularity  to  be  removable, 
the  following  limit  must  exist. 

lim  f{0) 

0-^0 
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Oi 


Insight  can  now  be  gained  by  looking  at  the  physical  nature  of  the  problem.  The 
singularities  occur  when  Tfe  =  0,  which  implies  that  e  =  0  in  Eq.  (3.6),  i.e.,  there  is  no 
uncertainty  in  that  measurement.  For  n  =  2,  this  would  mean  zj^  =  and  a  solution 
can  be  found  by  simple  least  squares,  i.e., 

lim/(«)=[HrHi]-'Hfzj  (3.10) 

Unfortunately,  this  formulation  is  only  valid  for  n  <  2,  because  for  n  >  2,  only 

has  rank  2  and  is  therefore  not  invertible. 

However,  one  can  expand  on  the  notion  of  a  noiseless  measurement,  and  use  the 
Kalman  filter  update  equations  to  add  an  additional  measurement  to  a  previous  estimate. 
For  n  =  2  and  n  =  3,  it  is  readily  apparent  that  only  one  T*,  can  be  equal  to  0  for  a  given 
experiment.  For  n  =  2,  the  singularities  are  points  corresponding  to  02  =  ~^k'>  since 
two  measurements  cannot  be  taken  at  the  same  frequency,  the  points  cannot  be  repeated. 
For  n  =  3,  the  individual  singularity  lines  lie  in  the  plane  02  =  — and  again,  since  no  two 
measurements  can  be  taken  at  the  same  frequency,  the  singularity  lines  cannot  intersect.  It 
is  suspected  that  this  is  the  case  for  higher  n  as  well. 
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This  implies  that  there  is  only  one  “noiseless  measurement”  at  any  given  point.  In 
other  words,  only  one  term  in  each  of  the  summations  is  undefined  at  any  given  singular 


point: 

These  terms  are  not  defined  because  =  0,  causing  =  0.  Taking  the  Kalman  filter¬ 
ing  approach,  let  the  pre-update  covariance  and  estimate  be  the  covariance  and  estimate 
calculated  by  using  all  the  terms  except  the  fcth  one. 

P  =  [ ^  HjR.-'Hfc]  and  0  =  P  ^ 

\k^k  / 

Then  following  the  formulation  in  Ref.  [14],  it  can  be  shown  that  in  any  neighborhood 
around  the  singularity, 

m  "I  r  m 

f{e)= 

.k=i  J  Lfe=i 

=  [p-'  +  Hl-Rj'Hj]"'  [p-'e  +  Hl-Rj'zi 

=  P-PH|’(Rj  +  HjPHf)''Hjp]  [p-'9  +  HfRj->2i 

=  pp-i0  -  PH?’  (Rfc  +  H^PH[)"^  H^PP-^0 

+PHrRr>.^  -  PH?  (Rj  +  HjPHr)^‘  HiPH^Rf  z; 

=  M-PHr(Ri  +  HtPH!')''HjW 

+pHr  [Rr>  -  (Rj  +  HiPHf)"'  HjpHfRf  ]  zj 

=  e-PH?(Rj  +  HjpHf)-‘Hi9 

+PHr  (Rj  +  HtPH?)-'  [(Rj  +  HiPH?)  Rf  -  HiPHJ-Rj']  zj 

=  9-pH?'(Rj  +  HiPHT)-'Hifl 

+PHr  (Rj  +  H,PhT)-‘  [RjRr>  +  HiPHTRj'  -  HjPHTRr']  zj 

=  e  -  PH?  (Rj  +  HiPHf)  ■'  Hje  +  PHl-  (Rj  +  HtPH?)  '■  Izj 

=  9  +  PHf  (Ri  +  HjPhT)"'  (zj  -  Hj9) 
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Therefore,  define  g{6)  to  be 


9(9)  =  9  +  PH?  (h^PH?  +  R()  (zi  -  Hj9) 

As  can  be  seen,  g{0)  does  not  have  a  singularity  when  Rj^  =  0,  and  for  any  given  k 
singularity,  f{0)  =  g{9)  in  a  deleted  neighborhood  around  that  singularity.  But  g(0) 
is  defined  and  continuous  on  the  entire  neighborhood,  so  extending  f(9)  by  filling  in  its 
singularities  with  g(9)  results  in  a  continuous  function  that  is  defined  and  continuous  on 
the  whole  of  7^^”. 


3.3.2  Boundedness.  It  is  now  necessary  to  show  that  f(9)  is  a  bounded  function. 
To  show  boundedness  for  a  continuous  function,  one  can  show  that  the  limit  as  9  grows 
unboundedly  exists  and  is  bounded.  First,  define 

n  =  {h£Tl'^^  :  II  h  11=  l} 

0  =  ah  where  h  €  and  a  G  7^  (3-11) 
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which  leads  to 


sup  II  f{9)  ||:=  sup  II  m(a)  ||=  5"  <  oo 

6  « 

Since  the  supremum  of  the  infinity  norm  of  /  is  5,  this  implies  that  /  maps  the  convex, 
bounded  set  5(0, 52")  c  7^^”  into  itself.  Now,  since  f{9)  is  a  continuous  function  that  maps 
a  convex  compact  set  into  a  convex  compact  set,  then  as  a  consequence  of  the  Schauder 
Fixed  Point  theorems  [16],  there  exists  a  point  9  G  5(0,52")  such  that  f{9)  =  9,  i.e.,  a 
fixed  point  exists. 
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3-4  Determination  of 


The  calculation  of  the  R  matrix  is  primarily  dependent  on  the  Tfc  matrix,  but  the 
R^b  matrix  is  also  required.  If  the  measurement  noise  is  on  Ak  and  Bk  directly,  then 
this  is  not  a  problem.  However,  to  make  it  more  complicated,  direct  measurements  of  Ak 
and  Bk  are  generally  not  available.  Rather,  Fourier  analyzer  supplied  measurements  of 
magnitude  and  phase  angle  are  available.  Therefore  in  this  section,  three  different  scenarios 
are  discussed;  1)  Constant  strength,  uncorrelated  noise  on  the  real  and  imaginary  phasor 
components  directly;  2)  Constant  strength,  uncorrelated  noise  on  the  magnitude  and  phase 
angle  in  radians  as  supplied  by  Matlab;  and  3)  Constant  strength,  uncorrelated  noise  on 
the  magnitude  in  decibels  and  phase  angle  in  degrees  as  supplied  by  a  frequency  analyzer. 

3.4- 1  Noise  on  the  Phasor  Components  Ak  and  Bk-  This  is  the  most  unlikely 
scenario  of  the  three,  but  it  is  also  the  easiest  and  the  one  on  which  the  derivation  in 
Section  3.2  is  based.  Here,  one  assumes  direct  access  to  the  real  and  imaginary  parts  of  the 
frequency  response  data.  This  implies  that  the  measurement  noise  strengths  and  correlation 
are  known,  i.e.,  R^b  is  known.  A  plot  of  a  representative  noise  is  shown  in  Fig.  3.2. 


Figure  3.2.  Constant  Strength  Uncorrelated  Noise  Added  to  Ak  and  Bk 
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3.4-2  Noise  on  the  Matlab  Supplied  Mk  and  (f>k.  As  stated  previously,  direct  mea¬ 
surements  of  Ak  and  Bk  are  generally  not  available.  The  more  likely  scenario  is  that  one 
has  hardware  or  software  supplied  measurements  of  magnitude  (Mk)  and  phase  angle  (<^fc). 
The  following  equations  are  then  used  to  calculate  Ak  and  Bk- 

Ak  =  Mk  cos  Bk  =  Mk  sin  (t>k  (3.12) 

It  is  more  likely  that  the  measurement  noises  {\m<i>)  on  the  observables  Mk  (A/'(0,<t^)) 
and  ^k  (A/’(0,(T^))  are  uncorrelated,  and  this  would  lead  to  correlation  between  and 
vb^.  Representative  noises  for  this  scenario  are  shown  in  Fig.  3.3.  As  can  be  seen  here, 
another  complication  is  that  the  strength  of  the  noise  on  Ak  and  Bk  is  different  at  each 
frequency.  To  remedy  this,  one  can  assume  that  the  noises  are  small,  and  using  Eq.  (3.12), 


Frequency 


Frequency 


Figure  3.3.  Constant  Strength  Noise  on  Mk  and  (pk 


the  following  linearized  relation  is  derived. 


AAfc 

■ 

cos  ^k 

-Mk  sin  (pk 

AMk 

ABk 

✓ 

sincpk 

Mk  cos  pk 

Apk 

V 

Now  Kab  is  calculated  as  follows. 

R-ab  =  E  {vfcv^ }  =  TM<t,E  (3.13) 
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This  indicates  that  R^ib  will  be  different  for  each  measurement,  and  will  need  to  be  re¬ 
calculated  each  time  before  adding  to  the  block  diagonal  R  matrix.  In  practice,  this 
calculation  is  not  completely  accurate,  since  the  true  values  of  Mk  and  which  are  re¬ 
quired  for  are  not  available.  One  can  use  the  measured  values  in  their  place  and 

achieve  good  results  for  smaller  noise  values,  but  much  better  results  are  obtained  if  one 
uses  the  parameter  estimates  to  estimate  the  transformation  matrix  as  well. 
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3.5  Other  Methods  of  System  Identification 

As  mentioned  in  Chapter  II,  there  are  a  couple  common  methods  of  dealing  with 
correlated  system  noises.  The  first  method  is  the  Instrumental  Variable  method.  In  this 
method,  one  calculates  the  next  estimate  using  an  H  matrix  which  is  the  H  matrix  that 
would  have  resulted  from  a  system  consisting  of  the  estimated  parameters.  In  other  words, 
the  estimated  parameters  yield  an  estimated  system.  The  elements  of  the  H  matrix  are  then 
the  transformed  magnitude  and  phase  values  off  of  the  estimated  Bode  plot.  The  following 
equation  is  then  iterated  until  a  suitable  convergence  criterion  is  met: 

=  (h’’h)“^H^z  (3.16) 

The  second  common  method  for  dealing  with  correlated  noise  is  the  Generalized 
Least  Squares  method.  A  purported  example  of  this  is  contained  in  Ref.  [19].  The  derivation 
of  the  noise  model  is  essentially  the  same,  but  too  many  simplifying  assumptions  are  made 
about  the  nature  of  the  measurement  noise.  The  final  result  can  be  shown  to  be  identical 
to  a  Generalized  Minimum  Variance  algorithm  that  assumes  =  I. 


3.6  Second  Order  Example 

For  investigative  purposes,  a  general  second-order  dynamical  system,  which  is  repre¬ 
sentative  of  an  aircraft’s  elevator-to-pitch  rate  transfer  function,  is  considered.  The  actual 
values  are  obtained  from  [2],  and  the  transfer  function  is  given  by 

T(s)  =  ^  =  4.8^  +  1.44 

u{s)  -  ais  -  a2  5^  -)-  0.84s  +  1.44  ^  ' 

The  Bode  plot  for  this  transfer  function,  which  is  representative  of  the  pitch  dynamics 
of  an  aircraft,  and  which  is  used  for  inner-loop  flight  control  system  design,  is  shown  in 
Fig.  3.5.  Also  shown  are  the  forty  measurement  frequencies  used  in  each  of  the  following 
experiments. 
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Figure  3.5.  Second  Order  Bode  Plot 
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As  can  be  seen,  if  a\  =  and  if  va  and  vb  are  uncorrelated,  then  reduces 

to  o'^l2x2-  However,  one  stiU  cannot  use  Eq.  (2.2)  to  estimate  6  because  Rj,  /  Rj  when 
^  • 

3.6.1  Basis  for  Choice  of  Noise  Strengths.  During  the  course  of  this  example,  the 
three  scenarios  discussed  in  Section  3.4  are  examined.  An  effort  is  made  to  keep  the  noise 
comparable  in  each  of  the  sections.  The  noise  strengths  used  in  this  example  are  representa¬ 
tive  of  a  Tektronix  frequency  analyzer.  The  manufacturer’s  specifications  give  measurement 
error  values  of  ±0.2  dB  and  ±0.5  deg.  These  are  taken  conservatively  as  two  sigma  values 
for  the  noise  on  the  amplitude  and  phase  measurements. 

For  the  constant  dB  and  deg  addition  in  scenario  2,  a  noise  with  a  covariance  matrix 
as  shown  in  Eq.  (3.15)  with  =  (0.2/2)^  and  =  (0.5/2)^  is  generated  and  added  to 
the  true  dB  and  degree  measurements  obtained  from  the  system  shown  in  Fig.  3.5.  These 
noisy  dB  magnitude  and  phase  values  are  then  converted  to  noisy  amplitude  and  phase 
values,  and  finally  to  noise- corrupted  A  and  B  values.  Each  of  the  noisy  values  is  then 
subtracted  from  its  true  value  to  obtain  the  transformed  noise.  The  results  of  this  addition 
and  transformation  process  are  shown  in  Fig.  3.4. 

In  Section  3.6.2,  constant  strength  noise  additions  to  Ak  and  Bk  are  examined.  To 
determine  a  reasonable  value  for  the  noise  strength  in  this  experiment,  the  actual  covariance 
of  the  Ak  and  Bk  noise  in  Fig.  3.4  is  calculated.  The  resulting  covariance  matrix  is  not  a 
scaled  version  of  the  identity  matrix,  but  for  this  experiment,  an  average  of  the  diagonal 
terms  is  used  for  and  ag,  and  the  off-diagonal  terms  are  set  to  zero.  This  results  in 
noise  as  shown  in  Fig.  3.2  with  a  covariance  matrix  of 

0.0232  0 

Kab  ~  (3.19) 

0  0.0232 

Section  3.6.3  discusses  the  identification  results  as  the  measurement  noise  strength 
is  increased  by  powers  of  10,  i.e., 

RABne^  =  l0^*'B.AB  p=  0,1, 2, 3, 4  (3.20) 
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This  higher  noise  evaluation  is  also  performed  for  the  other  noise  scenarios  discussed  in 
Section  3.4. 

In  Section  3.6.5,  a  constant  strength  noise  is  added  to  the  amplitude  and  phase 
variables.  To  determine  appropriate  strengths  for  these  noises,  the  covariances  of  the  Mk 
and  (f)k  noises  shown  in  Fig.  3.4  are  calculated,  and  these  covariances  are  used  to  form 
constant  strength  noise  additions  as  shown  in  Fig.  3.3.  The  resulting  covariance  matrix  is 

0.03172  Q 

=  (3.21) 

0  0.00442 

3.6.2  Inadequacy  of  the  Least  Squares  Method.  In  this  section,  one  pretends  that 
Ak  and  Bk  are  directly  measurable.  To  compare  the  different  estimation  methods,  uncorre¬ 
lated  and  equal  strength  noises  are  added  to  the  true  Ak  and  Bk  after  they  are  computed 
from  a  Bode  analysis  of  the  transfer  function  in  Eq.  (3.17).  Representative  noise  is  shown 
in  Fig.  3.2  with  the  resulting  values  of  =  ag  =  0.023.  This  choice  of  results  in  a 
diagonal  R  matrix,  but  it  is  stiU  not  a  scalar  multiple  of  the  identity  matrix  because  of  the 
varying  measurement  frequencies. 

Figure  3.6  displays  the  results  of  a  100  run  Monte-Carlo  (MC)  analysis,  which  upon 
evoking  the  weak  law  of  large  numbers  [11],  renders  a  gauge  of  the  identification  algorithm’s 
estimation  bias.  The  estimation  results  are  first  normalized  by  dividing  each  estimate  by  the 
true  estimate,  and  then  plotted.  Ellipses  are  plotted  representing  each  estimation  method’s 
actual  one  sigma  variation.  The  ellipses’  axes  intersect  at  the  average  estimate  for  each 
method.  The  algorithm  predicted  estimation  error  covariances  are  not  plotted,  because 
they  are  not  known  for  the  LS  or  GLS  methods,  i.e.,  Eq.  (2.2)  cannot  be  calculated  because 
there  is  no  constant  a  that  can  be  used.  Indeed,  the  cr  required  in  Eq.  (2.2)  is  estimated  as 
follows: 
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Figure  3.6.  All  Estimates  with  Uncorrelated  Ak  and  Bk  Noise 

However,  this  estimate  is  data  driven  and  the  predicted  estimation  error  covariance  stiU 
would  not  be  an  accurate  representation  of  the  true  error  covariance.  The  available  algo¬ 
rithm  predicted  covariances,  along  with  all  other  numerical  results,  are  given  in  Table  B.l. 

As  can  be  seen  in  the  plots,  there  are  large  biases  in  the  LS  average  estimate,  and 
the  true  parameter  is  well  outside  the  one  sigma  bounds.  This  indicates,  as  is  shown  in  the 
plot,  that  the  majority  of  the  parameter  estimates  in  the  hundred  runs  are  further  from  the 
true  estimate  than  one  sigma. 

Dramatic  improvements  can  be  made  by  using  one  of  the  other  estimation  methods. 
As  seen  in  the  plot,  the  IV  method  does  have  quite  a  large  sigma  value,  but  the  estimate 
is  not  nearly  so  biased  as  the  LS.  The  GLS  and  GMV  estimates  cannot  be  distinguished 
on  this  plot  because  they  are  tightly  clustered  around  the  true  estimate.  Figure  3.7  zooms 
in  to  give  a  better  look  at  these  estimates.  This  magnification  provides  a  good  view  of  the 
IV  estimates  and  sigma,  but  the  GLS  and  GMV  estimates  are  stiU  clustered  rather  tightly 
around  the  true  parameter.  Indeed,  this  is  always  the  case.  The  GLS  and  GMV  estimates 
are  always  better  than  the  LS  and  IV.  For  this  reason,  the  LS  and  IV  estimates  are  dropped 
from  most  plots  in  the  sequel.  The  numerical  results  are  stiU  included  in  the  appendix. 
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Figure  3.7.  Magnification  of  Estimates  with  Uncorrelated  Ak  and  Noise 


Further  magnification  in  this  noise  scenario  would  provide  little  information  on  the 
comparison  of  the  GLS  and  GMV  estimates  because  they  are  identical  in  this  case.  The 
only  difference  between  the  two  is  a  scaling  of  the  R  matrix  which  cancels  out  in  the  final 
equation.  The  only  thing  gained  by  the  GMV  method  is  an  accurate  estimation  of  the 
estimation  error  covariance.  As  is  seen  in  Table  B.l,  these  covariance  estimates  are  rather 
close  to  the  actual  ones.  The  small  biases  and  covariance  differences  imply  that  the  GMV 
and  GLS  system  identification  algorithms  yield  relatively  unbiased  parameter  estimates. 
Additionally,  GMV  is  doing  a  good  job  of  predicting  the  accuracy  of  its  estimate,  I'.e.,  the 
algorithm  is  “efficient.” 

3.6.3  Increasing  Noise.  The  next  portion  of  the  experiment  involves  increasing  the 
level  of  the  measurement  noise,  as  in  Eq.  (3.20).  The  results  for  p  =  1  are  very  similar  to 
p  =  0.  But  when  p  is  increased  to  2,  a  trend  that  causes  problems  later  on  starts  to  appear. 
As  the  strength  of  the  measurement  noise  increases,  the  LS  estimates  start  to  cluster  around 
the  origin  (Fig.  3.8).  This  may  seem  like  it  should  not  effect  the  GMV  estimate,  and  it 
does  not  in  this  case,  but  recall  that  the  GMV  estimate  is  initialized  with  this  LS  estimate. 
When  p  is  increased  to  3,  the  GMV  estimates  also  start  to  cluster  around  0.  This  is  not  a 
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Figure  3.8.  Higher  Noise  Estimates  with  Uncorrelated  Ak  and  Bk  Noise 


gradual  migration  like  the  one  that  happens  in  the  LS  estimates  as  the  noise  increases;  it  is 
sudden.  Indeed,  it  is  caused  by  the  appearance  of  a  second  fixed  point  close  to  the  origin. 

To  better  illustrate  this  problem,  a  “shotgun  approach”  is  taken  for  one  particular 
noise  realization.  The  results  of  the  GMV  algorithm,  when  it  is  initialized  with  a  series  of 
different  vectors  whose  elements  range  from  -0.5  to  0.5,  are  shown  in  Figure  3.9.  The  left 
figure  shows  the  paths  taken  during  the  iterations.  The  first  iteration  step  is  shown  by  a 
dotted  line,  while  the  remaining  iterations  are  plotted  as  a  solid  line.  The  end  points  (or 
fixed  points)  are  shown  as  an  *.  Most  of  the  initial  value  choices  converge  to  the  desired 
estimate,  but  several  runs  initialized  close  to  the  origin  converge  to  a  different  value.  This 
seems  to  be  a  detrimental  property  of  the  phasor  oriented  GMV  algorithm  at  high  noise. 
How  does  one  choose  the  “correct”  fixed  point?  The  solution  to  this  may  lie  in  the  right  hand 
plot.  This  is  a  plot  of  [f{0)  —  O]^[f{0)  —  0]  at  each  iteration.  The  slowly  converging  paths 
in  this  plot  correspond  to  those  paths  that  are  converging  to  the  wrong  point.  Although  the 
LS  estimate  is  a  convenient  point  to  start  the  GMV  algorithm,  it  is  not  the  wisest  choice, 
since  it  tends  to  approach  the  origin  as  the  noise  level  increases. 
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Figure  3.9.  Convergence  Properties  for  Uncorrelated  Ak  and  Bk  Noise 


A  better  choice  would  be  a  point  in  the  same  quadrant  as  the  true  parameter,  but 
with  a  larger  magnitude.  This  is  usually  not  possible,  because  one  does  not  know  the  true 
parameter.  However,  if  the  problem  is  known  to  be  bothered  by  high  noise,  one  can  choose 
several  starting  points,  perhaps  one  from  each  quadrant,  and  the  one  that  converges  the 
fastest  is  the  correct  one. 

If  one  examines  the  convergence  rates  for  the  p  =  3  estimates  in  Fig.  3.8,  which  are 
shown  in  Fig.  3.10,  it  can  be  seen  that  many  estimates  converge  very  slowly,  while  others 
do  not  appear  to  converge  at  aU.  However,  when  the  GMV  algorithm  is  initialized  at  the 
point  00  =  [—100  —  100  0  0]^,  the  results  improve  dramatically.  These  results  are  shown 
in  Fig.  3.11.  As  can  be  seen,  the  GMV  estimates  are  no  longer  trapped  around  the  origin, 
and  the  convergence  rates  are  much  quicker  as  well. 

3.6.4  Noise  on  Observables.  In  this  section,  the  noise  is  added  to  the  dB  magnitude 
and  phase  angle  (in  degrees),  as  is  given  in  the  Tektronix  [33]  specifications.  Then,  the  two 
required  transformations  (Eqs.  (3.13)  and  (3.15))  are  performed  to  obtain  the  correct 
(Tg,  and  <r\g  used  for  the  matrix  calculation  in  Eq.  (3.18).  Through  proper  modeling 
of  the  noise  transformations,  the  GMV  algorithm  achieves  better  estimates  than  any  of  the 
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Iteration 


Figure  3.10.  Convergence  Rates  for  LS  Initialization 


Figure  3.11.  Estimates  and  Convergence  Rates  for  Modified  Initialization 


other  methods.  The  GLS  algorithm  comes  the  closest,  and  a  comparison  of  these  for  the 
p  =  0  case  is  shown  in  Fig.  3.12. 

p=0  p=0 


Figure  3.12.  Estimates  with  Uncorrelated  and  (t>deg  Noise 

As  can  be  seen,  the  GMV  estimates  have  both  a  lower  bias  and  error  covariance  than 
the  GLS  estimates.  The  reason  for  this  is  the  proper  modeling  of  the  matrix.  Initial 
attempts  at  this  used  the  measured  values  in  the  transformation  matrices.  While  this  yields 
good  results  for  small  noise  strengths,  the  measured  values  for  higher  noise  strengths  do  not 
provide  an  accurate  estimate  of  the  magnitude  and  phase  required  for  the  transformation 
matrix  calculations.  Problems  arise  in  poor  numerator  estimates,  and  the  parameters  start 
to  get  trapped  around  the  origin  sooner  than  they  would  otherwise. 

Instead,  the  measured  values  are  used  only  for  the  first  five  GMV  iterations.  This  is 
only  done  to  get  a  more  accurate  estimate  of  the  Bode  plot,  and  may  not  even  be  necessary. 
After  the  fifth  iteration,  the  GMV  parameter  estimates  are  used  to  estimate  the  Bode  plots, 
and  the  magnitude  and  phase  values  from  this  estimation  are  used  in  the  transformation 
matrices.  This  method  produces  much  more  accurate  estimates,  and  it  is  also  more  robust 
to  the  ’’zero  trapping”  than  before. 
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Higher  noise  strength  runs  were  also  performed  for  this  scenario,  and  the  numerical 
results  are  given  in  Table  B.3.  As  before,  the  GMV  results  are  the  best,  and  the  error 
covariance  prediction  is  close  to  the  actual. 

3.6.5  Noise  on  Amplitude  and  Phase.  In  this  section,  the  noise  on  the  amplitude 
and  phase  is  assumed  to  be  of  constant  strength  and  temporally  uncorrelated.  To  this  end, 
the  covariance  matrix  in  Eq.  (3.21)  is  used  to  generate  noises  for  addition  to  Mk  and  0^. 
This  implies  that  only  one  transformation  needs  to  be  performed  to  determine  (t\,  a g,  and 
(T\g.  Representative  noises  for  this  section  are  shown  in  Fig.  3.3.  Once  again,  the  measured 
values  are  used  in  Tjv/^  for  the  first  five  GMV  iterations,  and  then  the  parameter  estimates 
are  used  to  estimate  Tm<^-  The  numerical  results  for  all  the  noise  strength  cases  are  given 
in  Table  B.2. 

One  concern  of  note  in  this  scenario  is  shown  in  Fig.  3.13.  This  is  a  plot  of  the  error 
mean  and  covariance  for  the  GMV  algorithm  as  each  measurement  is  added.  The  algorithm 
estimated  covariance  is  close  to  the  actual  covariance  after  about  22  measurements,  but  the 
actual  covariance  appears  to  begin  to  diverge  toward  the  end.  The  22  measurements  is  not 
of  concern,  because  this  is  highly  dependent  on  the  order  in  which  the  measurements  are 
added  to  the  linear  regression.  However,  the  cause  of  the  divergence  is  not  known.  It  has 
something  to  do  with  the  transformation,  because  the  covariance  increase  does  not 

appear  in  the  GLS  estimates.  A  possible  cause  is  the  negative  measured  magnitudes  at  this 
noise  level  in  this  scenario.  This  may  violate  some  physical  property  or  assumption,  but  a 
definitive  cause  can  not  be  determined. 

3. 7  Discussion 

The  normalized  numerical  results  of  all  the  experiments  are  summarized  in  Ta¬ 
bles  B.l  through  B.3  in  Appendix  B.  In  all  cases,  the  values  shown  correspond  to  a  forty 
measurement  linear  regression.  The  Monte-Carlo  averaged  estimation  error  (e)  and  sigma 
((Tg)  are  given  for  aU  cases,  and  the  algorithm  predicted  sigma  ((Tp)  is  given  for  the  MV  and 
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Figure  3.13.  Estimation  Performance  for  Increasing  Measurements 


IV  experiments.  In  each  case,  the  bias  in  the  MV  estimate  is  about  two  orders  of  magnitude 
smaller  than  in  the  naive  LS  estimate,  and  the  estimation  error  covariance  is  much  smaller 
as  well. 

For  p  <  2,  aU  100  of  the  estimated  Bode  plots  are  virtually  indistinguishable.  Ap¬ 
pendix  B  contains  the  final  estimation  results  for  p  =  2, 3.  The  top  plots  show  the  true  Bode 
plot  and  the  average  of  the  estimated  Bode  plots,  along  with  the  range  of  the  estimated 
Bode  plots.  As  can  be  seen,  the  spread  is  quite  small  even  for  p  =  2.  It  appears,  however, 
that  the  admittedly  small,  identification  provided,  uncertainty  is  well  within  the  system’s 
bandwidth,  i.e.,  it  is  structured  uncertainty. 

The  bottom  plots  show  the  true  poles  and  zero,  along  with  the  estimated  poles  and 
zeros.  For  p  <  2,  the  estimated  poles  lie  practically  on  a  horizontal  line,  and  even  at  p  =  2, 
the  vertical  variation  is  very  small,  i.e.,  the  identification  algorithm  renders  a  very  accurate 
estimate  of  the  damped  natural  frequency.  A  noticeable  bias  begins  to  appear  when  p 
reaches  3,  but  it  is  not  all  that  large,  considering  the  noise  level. 
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3. 8  Conclusions 


In  this  chapter,  a  frequency  domain  approach  is  taken,  and  phasors  are  used  for 
system  identification.  Gaussian  measurement  noise  is  assumed,  as  is  customary  in  classical 
filtering  and  system  identification  work.  The  proper  minimum  variance  estimate  equations 
are  derived  and  applied  to  a  nonlinear  estimation  problem.  The  results  are  then  compared 
to  the  simple  minded  least  squares  estimate  for  a  second-order  system  that  is  representative 
of  an  aircraft’s  pitch  dynamics,  which  is  used  for  inner-loop  flight  control  system  design. 

The  GMV  estimate  outperformed  the  LS  estimate  in  all  cases.  The  LS  estimate  did, 
however,  provide  a  useful,  albeit  dangerous,  starting  point  for  iterating  the  GMV  estimate. 
At  small  noise  levels,  it  does  not  matter  where  the  GMV  algorithm  is  initialized,  but  as 
higher  noise  levels  appear,  it  is  dangerous  to  initialize  the  algorithm  close  to  the  0  point. 
There  is  an  additional  fixed  point  appearance  there  that  has  proven  capable  of  trapping  the 
estimate.  A  positive  note,  when  the  estimate  is  trapped,  the  algorithm  takes  much  longer 
to  converge,  so  multiple  high  magnitude  starting  points  should  allow  one  to  find  the  desired 
fixed  point.  In  conclusion,  the  least  squares  estimate  is  not  an  effective  one,  even  in  cases 
of  small  measurement  noise.  Although  the  IV  estimate  is  somewhat  better  than  the  LS 
estimate,  it  still  performs  poorly  under  increased  noise  levels.  The  GLS  and  GMV  provide 
much  more  accurate  estimates,  and  the  GMV  seems  to  provide  lower  biased  and  smaller 
covariance  estimates  in  each  case. 

In  this  chapter,  it  is  shown  that  the  frequency  domain  can  be  used  in  system  identi¬ 
fication.  Careful  stochastic  modeling  of  the  estimation  problem  at  hand  renders  an  efficient 
identification  algorithm  that  is  superior  to  straightforward  least  squares.  However,  the 
experimental  phasor  approach  is  applicable  to  stable  systems  only. 
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IV.  Effects  of  Unmodeled  Dynamics  on  the  Phasor  Approach 


When  system  identification  is  performed  for  control  system  design,  it  is  oftentimes 
required  to  identify  the  plant’s  dominant  mode  only.  Thus,  identification  is  performed 
in  the  presence  of  additive  modes,  and  one  then  refers  to  unmodeled  dynamics.  In  this 
chapter,  the  proper  choice  of  identification  signals  when  unmodeled  dynamics  are  present 
is  investigated.  Different  combinations  of  high  and  low  frequency  dynamics  are  added 
to  a  representative  and  physically  motivated  second-order  plant,  and  carefully  controlled 
identification  experiments  are  performed.  Moreover,  in  adding  additional  dynamics,  one  is 
strongly  guided  by  physical  considerations,  and  this  results  in  a  well-designed  experiment. 
To  solidify  the  concept,  consider  the  flight  control  context. 

Flight  control  system  designers  often  use  a  plant  model  which  consists  of  the  short 
period  pitch  dynamics  of  an  aircraft.  However,  the  short  period  dynamics  are  not  the  only 
modes  present  in  the  aircraft’s  pitch  channel.  This  directly  impacts  system  identification. 
Thus,  it  is  desired  to  identify  the  dominant  mode  of  a  system  when  there  are  other  un¬ 
modeled  dynamics  present  as  well.  Now,  the  more  one  knows  about  a  particular  system  a 
priori,  the  easier  the  identification  problem  becomes.  Prior  information,  including  band¬ 
width  information,  is  at  a  premium  in  system  ID.  Hence,  the  focus  in  this  chapter  is  on 
the  interaction  between  the  unmodeled  dynamics  and  the  proper  choice  of  input  frequencies 
used  to  identify  a  plant. 

In  Section  4.1,  a  simple  second-order  plant  representative  of  the  short  period  dynam¬ 
ics  of  an  aircraft,  and  no  additional  dynamics,  is  used  to  establish  a  baseline  for  determining 
optimal  excitation  frequencies.  The  short  period  dynamics  are  routinely  used  in  flight  con¬ 
trol  design  work.  Low  and  high  frequency  modes  are  added  to  the  second-order  plant  in 
Sections  4.2  and  4.3,  respectively.  The  idea  here  is  that  modeling  based  on  physical  con¬ 
siderations  renders  mathematically  tractable  problems.  Hence,  the  low  frequency  mode  is 
chosen  to  represent  the  aircraft’s  phugoid  mode,  and  the  high  frequency  dynamics  model 
a  flexible  mode.  Both  the  low  and  high  frequency  modes  are  included  in  the  analysis  and 
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experiments  presented  in  Section  4.4,  and  the  restdts  of  naively  overmodeling  a  system  are 
in  Section  4.5. 


4.1  Second-Order  Plant 


In  this  section,  the  proper  choice  of  frequencies  used  to  identify  a  second-order 
underdamped  plant,  with  one  zero,  is  examined.  This  constitutes  a  baseline  for  comparing 
with  results  obtained  in  the  presence  of  unmodeled  dynamics.  For  this  purpose,  the  plant 
given  in  Eq.  (4.1),  which  is  representative  of  the  short  period  dynamics  of  an  aircraft,  is 
used.  The  Bode  plot  of  its  transfer  function  is  shown  in  Fig.  4.1,  from  which  the  plant’s 
bandwidth  is  readily  apparent. 

4.8s  -1-  1.44  4.8(s  +  0.3) 


T2(s)  = 


s2  +  0.84s  +  1.44  s2  -I-  2(0.35)(1.2)s  +  (1.2)2 


(4.1) 


Figure  4.1.  Second-Order  Bode  Plot 


To  examine  the  impact  on  identification  of  the  choice  of  input  frequency,  the  transfer 
function  in  Eq.  (4.1)  is  used  to  generate  a  series  of  measurements,  and  Gaussian  noise  of 
strength  (TMaB  —  ^<t>deg  —  is  then  added  to  the  measurements.  The  frequency 
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range  shown  in  Fig.  4.1  is  then  divided  to  provide  16  measurements  per  half  decade.  Finally, 
a  half  decade  window  is  stepped  through  the  given  frequency  range  at  1/4  decade  incre¬ 
ments,  resulting  in  15  different  parameter  estimation  frequency  windows.  For  each  window, 
twenty  Monte- Carlo  runs  are  performed  using  the  estimation  algorithm  described  in  the 
previous  chapter.  These  results  are  graphically  shown  in  Fig.  4.2,  and  the  actual  values  are 
summarized  in  Table  C.l. 


al  Estimation  Error 


a2  Estimation  Error 


Measurement  Window 


Measurement  Window 


bl  Estimation  Error  b2  Estimation  Error 


Measurement  Window  Measurement  Window 


Figure  4.2.  Results  of  Identifying  the  Second-Order  Plant 


As  can  be  seen,  the  bias  in  the  estimate  is  extremely  large  when  very  low  frequencies 
are  used  for  identification.  An  interesting  note  is  that  the  covariance  is  small.  This  is 
truly  poor  ID  because  the  algorithm  believes  that  its  estimate  is  accurate,  when  in  fact  the 
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opposite  is  true.  In  the  parlance  of  Kalman  Filtering,  one  refers  to  filter  divergence  -  an 
unfortunate  manifestation  of  nonlinearities  encountered  when  extended  Kalman  filters  are 
used.  Furthermore,  the  algorithm  is  identifying  the  plant  as  a  first-order  unstable  plant, 
which  may  be  caused  by  the  effects  of  the  transfer  function  zero  at  lower  frequencies.  The 
covariance  then  increases  a  bit  as  the  bias  approaches  zero,  and  returns  to  small  values 
around  the  center  of  the  frequency  range. 

The  estimation  bias  is  smallest  overall  at  window  number  seven,  which  corresponds 
to  the  half  decade  window  just  below  1  rad/sec.  This  window  is  one  below  that  which 
contains  the  second-order  mode’s  natural  frequency.  Also,  the  pole-zero  estimation  is  most 
accurate  at  the  same  window.  The  estimation  bias  then  becomes  large  again  at  the  high 
frequencies,  but  the  covariance  increases  as  well.  This  is  a  better  result  than  that  yielded  by 
the  low  frequency  inputs  because  the  algorithm  is  warning  the  user  that  its  estimate  is  poor. 
The  sole  exception  to  the  increasing  bias  is  the  estimate  of  6i,  which  remains  accurate.  This 
seems  to  imply  that  the  gain  estimate  may  stiU  be  accurate. 

The  low  and  high  measurement  frequency  problems,  noted  in  this  section  for  a 
second-order  plant,  form  the  basis  for  the  analysis  in  the  following  sections.  These  sections 
include  the  addition  of  low  and/or  high  frequency  dynamics. 

4.2  Low  Frequency  Unmodeled  Dynamics 

In  this  section,  a  second-order  identification  is  again  performed.  However,  the  plant 
dynamics  here  include  an  unmodeled,  low  frequency,  mode  in  addition  to  the  short-period 
mode.  Guided  by  insights  into  the  physics  of  flight  dynamics,  this  mode  is  similar  to  the 
phugoid  mode  of  an  aircraft  [2].  Thus,  the  fourth-order  elevator  deflection  to  pitch  rate 
transfer  function  is  given  by 

T  (s)  =  4.7843s(a-b0.016)(a  + 0.2862) 

’  (s2 0.004665 -I- 0.0053)(s2 -f  0.8394s  +  1.4383)  ^  ^ 

To  derive  T4i{s),  a  particular  low  frequency  mode  is  chosen  and  included  in  a  fourth- 
order  transfer  function.  Then,  a  numerical  optimization  routine  is  used  to  minimize  the 
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distance  between  the  Bode  plot  of  the  fourth-order  transfer  function  and  the  Bode  plot  of 
T2{s),  at  frequencies  between  and  10®'®  rad/sec.  The  transfer  function  in  Eq.  (4.1) 

results.  This  match  is  performed  in  a  noiseless  environment,  so  the  parameters  in  Eq.  (4.1) 
are  used  to  assess  the  error  in  the  estimate  obtained  assuming  a  second-order  model  and 
using  the  fourth-order  plant’s  measurements.  The  Bode  plots  of  the  plants  in  Eqs.  (4.1) 
and  (4.2)  are  compared  in  Fig.  4.3.  As  is  seen,  the  Bode  plots  match  around  the  natural 
frequency  of  the  “dominant”  mode  and  above,  but  differ  dramatically  around  the  natural 
frequency  of  the  low  frequency  mode  and  below. 


Figure  4.3.  Second-  and  Fourth-Order  Bode  Plot  Comparison 

The  same  moving  frequency  window  estimation  is  performed  using  the  noise  cor¬ 
rupted  outputs  of  the  fourth-order  plant.  The  graphical  results  are  portrayed  in  Fig.  4.4 
and  the  actual  numerical  values  are  contained  in  Table  C.2. 

In  this  case,  the  low  frequency  estimation  error  covariances  have  collapsed  to  approx¬ 
imately  zero,  while  the  bias  remains  very  large.  The  addition  of  the  low  frequency  dynamics 
has  hurt  the  estimation  process.  If  one  is  naive  enough  to  take  low  frequency  measurements 
into  account,  the  algorithm  would  produce  a  seemingly  accurate  estimate,  when  in  fact  it 
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Error  Mean  and  Covariance  Error  Mean  and  Covariance 


al  Estimation  Error 


bl  Estimation  Error 


a2  Estimation  Error 


Measurement  Window 


b2  Estimation  Error 


Figure  4.4.  Results  of  Identifying  the  Fourth-Order  Low  Frequency  Plant 


is  totally  wrong.  The  error  covariances  then  balloon  as  the  estimation  window  moves  into 
the  range  of  the  desired  mode  and  collapse  once  again  as  the  bias  approaches  zero.  Addi¬ 
tionally,  the  high  frequency  results  are  virtually  identical  to  those  of  the  previous  section. 
An  interesting  note  is  that  the  most  accurate  estimate  of  the  parameters,  poles,  and  zero 
has  now  moved  to  window  number  eight;  the  one  containing  the  desired  dominant  mode’s 
natural  frequency. 

4.3  High  Frequency  Unmodeled  Dynamics 

In  this  section,  the  procedure  is  repeated  with  the  low  frequency  mode  replaced  by  a 
high  frequency  mode.  Again,  the  problem  is  formulated  by  deliberately  leaning  toward  re¬ 
liance  on  physical  insight  into  the  problem.  Thus,  this  mode  is  similar  to  the  high  frequency 
bending  mode  of  an  aircraft  [26] ,  and  the  fourth-order  transfer  function  is  given  by 

5.733(s  4-0.2979X52  + 0.4965 +  695) 

(52  +  0.83945 +  1.4383)(52  + 2.865 +  827)  ^  ^ 

Here  a  high  frequency  mode  is  chosen  and  included  in  a  fourth-order  transfer  func¬ 
tion.  Then  an  optimization  routine  is  once  again  used  to  minimize  the  distance  between  the 
Bode  plot  of  the  fourth-order  transfer  function  and  that  of  Eq.  (4.1)  at  frequencies  between 
10“°'®  and  10°'®  rad/sec.  Therefore  as  before,  the  parameters  in  Eq.  (4.1)  are  deemed  to  be 
the  desired  parameters.  The  Bode  plots  of  the  plants  in  Eqs.  (4.1)  and  (4.3)  are  compared  in 
Fig.  4.5.  As  is  seen,  the  Bode  plots  match  at  low  frequencies,  but  differ  in  the  neighborhood 
of  the  high  frequency  mode. 

The  same  moving  frequency  window  estimation  is  performed  using  the  noise  cor¬ 
rupted  outputs  of  the  fourth-order  plant.  The  graphical  results  are  portrayed  in  Fig.  4.6 
and  the  actual  numerical  values  are  contained  in  Table  C.3. 

Once  again,  the  additional  mode  introduces  unmodeled  dynamics,  which  potentially 
hamper  the  identification  process.  The  low  frequency  performance  is  similar  to  the  second- 
order  results,  but  the  high  frequency  error  covariances  collapse  again  in  the  area  around  the 
additional  mode’s  natural  frequency.  Even  the  estimate  of  the  61  parameter  now  has  a  large 


Figure  4.5.  Second-  and  Fourth-Order  Bode  Plot  Comparison 

bias  and  small  error  covariance.  The  area  of  covariance  coUapse  is  not  as  large  as  when  the 
low  frequency  mode  is  added,  but  as  in  the  case  of  low  frequency  unmodeled  dynamics,  it 
corresponds  to  the  area  where  the  second-order  Bode  plot  differs  from  the  fourth-order. 

The  window  of  best  estimation  is  also  more  ambiguous.  Best  parameter  and  pole- 
zero  estimates  range  between  windows  eight  and  nine,  depending  on  the  parameter.  This 
seems  a  bit  strange  when  one  considers  the  high  frequency  mode  should  push  the  better 
identification  to  lower  frequencies.  However,  the  fourth-order  transfer  function  is  matched 
to  a  second-order  in  a  subset  of  window  eight,  which  implies  the  best  match  should  always 
be  in  window  eight. 

4.4  High  and  Low  Frequency  Unmodeled  Dynamics 

In  this  section,  both  the  high  frequency  and  low  frequency  modes  are  added  to  the 
plant,  and  the  system  identification  procedure  is  repeated.  This  is  similar  to  the  situation 
in  a  real  aircraft,  and  the  sixth-order  transfer  function  is  given  by 

,  _  5.7141s(s  -I-  0.016)(s  -b  0.2841)(s2  +  0.4965  +  695)  . 

~  (52  +  0.004665  -f  0.0053)(52  -f-  0.84245  -b  1.4405)(52  -b  2.865  -b  827)  ^  ’ 


4-8 


Error  Mean  and  Covariance  Error  Mean  and  Covariance 


al  Estimation  Error 


bl  Estimation  Error 


Figure  4.6.  Results  of  Id 
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Te{s)  is  the  result  of  using  the  same  high  and  low  frequency  modes  as  before,  and 
an  optimization  routine  to  match  a  sixth-order  plant’s  Bode  plot  to  that  of  Eq.  (4.1)  at 
frequencies  between  and  10°  ®  rad/sec.  Therefore  as  before,  the  parameters  of  this 

second-order  transfer  function  are  used  as  the  desired  parameters.  The  Bode  plots  of  the 
plants  in  Eqs.  (4.1)  and  (4.4)  are  compared  in  Fig.  4.7.  Obviously,  the  Bode  plots  now  differ 
at  both  low  and  high  frequencies. 


Figure  4.7.  Second-  and  Sixth-Order  Bode  Plot  Comparison 

The  same  moving  frequency  window  estimation  is  performed  using  the  noise  cor¬ 
rupted  outputs  of  the  sixth-order  plant.  The  graphical  results  are  portrayed  in  Fig.  4.8  and 
the  actual  numerical  values  are  contained  in  Table  C.4. 

The  estimation  problems  noted  before  are  now  compounded.  The  error  covariance 
is  now  artificially  small  at  both  low  and  high  frequencies.  This  collapse  of  estimation 
performance  occurs  when  the  Bode  plots  of  the  actual  and  desired  plants  differ  greatly. 
Additionally,  the  best  estimation  is  now  definitely  achieved  in  window  eight.  This  further 
supports  the  notion  that  the  best  estimate  occurs  in  the  window  where  the  Bode  plots  are 
matched. 
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Figure  4.8.  Results  of  Identi 
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the  Sixth- Order  Plant 


4-5  Attempt  at  Overmodeling 


In  this  section,  the  effect  of  overmodeling  in  the  face  of  limited  measurements  is 
examined.  The  plant  is  modeled  as  sixth-order,  but  each  estimate  stiU  only  uses  16  mea¬ 
surements.  As  expected,  the  results  are  very  poor,  and  representative  samples  are  shown 
in  Fig.  4.9.  Compounding  the  problem  that  16  measurements  are  not  sufficient  to  identify 
10  parameters  accurately,  is  the  fact  that  none  of  the  frequency  windows  provides  enough 
information  about  the  full  bandwidth  for  an  accurate  estimate.  This  emphasizes  the  idea 
of  keeping  the  number  of  estimated  parameters  small  as  related  to  the  measurement  infor¬ 
mation. 


al  Estimation  Error  a2  Estimation  Error 


Figure  4.9.  Overmodeling  Results 
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4-6  Conclusions 


Some  degree  of  a  priori  knowledge  about  a  plant  must  exist  for  accurate  identifica¬ 
tion.  In  the  context  of  this  work,  this  entails  knowledge  of  the  plant’s  bandwidth,  and  of 
the  bandwidth  of  the  dominant  mode. 

When  one  attempts  to  obtain  a  second-order  model  for  a  second-order  plant,  the 
measurements  should  be  in  the  region  of  the  mode  present  in  the  plant  or  an  accurate 
estimate  wiU  not  be  possible.  In  the  discussed  flight  control  scenario,  it  seems  that  the 
optimum  frequency  range  for  proper  estimation  is  immediately  prior  to  the  plant’s  natural 
frequency.  However,  this  may  be  specific  to  this  plant  and  its  low  frequency  zero. 

When  one  attempts  to  obtain  a  second-order  model  of  a  higher  order  plant,  i.e.,  when 
identification  in  the  presence  of  modeling  error  is  performed,  the  importance  of  applying 
sinusoidal  test  signals  whose  frequency  is  near  the  frequency  of  the  desired  mode  becomes 
greater.  Otherwise,  not  only  will  the  estimation  process  give  you  an  inaccurate  (biased) 
estimate,  but,  in  addition,  it  wiU  tell  you  that  it  is  an  accurate  estimate.  This  would  lead 
to  placing  undue  confidence  in  the  estimate  and  is  akin  to  the  notorious  filter  divergence 
phenomenon  encountered  in  extended  Kalman  filtering. 

The  best  parameter  match  seems  to  be  in  window  number  eight  in  all  higher  order 
cases.  This  is  probably  due  to  matching  the  second-order  Bode  plot  to  the  higher  order 
one  in  that  frequency  range.  If  a  different  range  had  been  used,  the  results  would  probably 
have  matched  in  that  range  instead.  Small  variances  in  the  matching  range  produce  small 
changes  in  the  parameters.  Therefore,  the  small  biases  in  the  windows  around  eight  are 
probably  due  to  matching  error  and  not  actual  estimation  error. 
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V.  Identification  of  a  Discrete-Time  Dynamical  System 


While  the  identification  of  continuous-time  plants  from  frequency  data  is  very  useful, 
the  primary  scenario  for  the  identification  of  dynamical  systems  lies  in  the  time  domain.  In 
the  control  systems  paradigm,  the  parameters  of  the  plant  are  needed  to  design  compen¬ 
sators  to  elicit  desired  responses  from  the  plant.  To  this  end,  plant  input  is  usually  applied 
in  the  present  day  by  computer  controller  actuators.  This  brings  the  discrete-time  plants 
models  into  play.  Control  is  typically  applied  through  a  zero-order  hold  device  to  the  plant. 

In  this  chapter,  the  discrete-time  parameters  of  a  given  physical  plant  are  identified 
using  several  methods.  A  comparison  is  made  between  the  proposed  Generalized  Minimum 
Variance  (GMV)  and  established  Generalized  Least  Squares  (GLS)  identification  algorithms. 
The  Instrumental  Variable  (IV)  method  is  not  presented  here,  because  the  results  are  similar 
to  those  observed  in  Chapter  III,  where  the  IV  estimates  were  significantly  worse  than  either 
the  GLS  or  GMV  estimates. 

Here,  the  experimental  data  consists  of  a  finite  number  of  input  and  output  mea¬ 
surements  of  a  discrete-time  plant,  which  is  representative  of  the  sampled  pitch  dynamics  of 
a  transport  aircraft.  Also,  sensor  noise  has  corrupted  the  measurements  of  the  output.  The 
GLS  method  is  described  in  Ref.  [14],  and  is  basically  the  same  as  the  method  proposed  in 
Ref.  [32]. 

Section  5.1  sets  up  the  time  domain  problem,  and  explains  the  application  of  the 
GMV  and  GLS  methods.  Section  5.2  presents  the  actual  plant  used  in  the  experiments.  The 
algorithms  are  compared  using  only  noise  on  the  output  in  Section  5.3.  The  convergence 
properties  are  discussed  in  Section  5.4.  Input  noise  is  added  to  the  problem  in  Section  5.5. 
Section  5.6  contains  the  results  when  the  sampling  rate  is  increased  by  a  factor  of  four. 
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5.1  Time  Domain  Identification  Setup 


Consider  the  discrete-time  2:-domain  transfer  function  given  by 

fr  \  ^  -f  b2Z^  ^  -h - \-  bn-iz  -|-  bn  _  J27=o  ^  /-g 

u(z)  2"  -  -  a22"-2  - - an-iZ  -  an 

where  ao  =  —1,  and  the  2n  -t-  1  parameters  ai . .  and  bo  -  ■  -  bn  are  unknown  and  may 
have  a  zero  value,  as  discussed  in  Chapter  III.  This  transfer  function  can  be  rewritten  in 
time  domain  form  as 

(q~^  —  q~^ai  — - vikT)  =  (q~^bo  -|-  q~^b2  + - q  '^bn'^  u{kT)  (5.2) 

^  **  —  1..^ 

where  T  is  the  sampling  period  and  q~‘^  is  an  n  step  time  delay.  Letting  pk  =  y{kT)  and 
Uk  =  u(kT),  Eq.  (5.2)  can  be  written  in  recursive  form  as  follows: 

Vk  =  O‘iyk-1  +  0.2yk-2  + - h  dnyk-n  +  +  biUk-1  +  b2Uk-2  + - 1"  Kuk-n  (5.3) 


To  identify  the  parameters  of  the  discrete  time  system,  one  need  only  solve  the  following 
system  of  2n  -f  1  equations 
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yk+3 
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2/fc-i 

yk—n-\-l 
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■  Uk—n+1 
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yk-2+N  •  • 
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Uk+N  ■ 
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bn 

However,  the  values  for  yk  and  Uk  cannot  be  directly  measured.  Rather,  we  can 
measure  the  noise  corrupted  Yk  and  Uk,  where  the  noise  is  assumed  to  be  Gaussian: 


Yk  =  yk  +  n,  V  =  Af{0,  al) 
Uk  —  Uk  d"  ^  ^^(0)  ^u;) 
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This  changes  Eq.  (5.3)  to 

Yk  Vh:  —  Cli(Yjf:—\  '^k—1^  "h  f^2(Yk—2  2)  "h  ■  '  ■  T  0,n(Yk—n  '^k—n) 

T  ^oiUk  ~  '^k)  T  1  1)  T  ■  ■  ■  T  bji(^Uk—n  ~  '^k—n)  (^•^) 

or 


Yk  —  ttllfc-l  +  a2Yk-2  +  •  •  •  +  dnYk-n  +  boUk  +  b\Uk-\  +  b2Uk-2  +  •  •  •  +  bnUk-n  +  Vk 


where 


Vk  =  Vk-  aiVk-1  -  a2Vk-2 - dnVk-n  “  boWk  -  biWk-1  -  b2Wk-2 - bnWk-n  (5.6) 


This  expression  is  now  in  a  form  that  can  be  set  up  in  a  statistical  Linear  Regression 
equation  as  in  Eq.  (2.7),  where 
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'  ■  Yk—n+l 

Uk+i 

Uk  ■ 

Uk-n-kl 
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[  Vk+N  \ 

As  before,  the  R  matrix  is  the  expected  value  of  ee^,  and  if  the  noise  is  assumed  to  be 
white,  then  the  R  matrix  will  be  an  A  x  A  Toeplitz  matrix  with  a  non-zero  diagonal  and 
n  non-zero  off  diagonal  terms  above  and  below  the  diagonal. 
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5.1.1  GLS  Estimation.  The  GLS  algorithm  is  derived  through  an  attempt  to  match 
the  identification  problem  to  a  rigorous  least  squares  formulation.  The  LS  estimate  would 
be  accurate  if  the  R  matrix  was  equal  to  a  scalar  multiple  of  the  identity  matrix.  This 
would  occur  if  the  series  was  white. 

Assuming  there  is  measurement  noise  only  on  the  output,  the  system  under  consid¬ 
eration  can  be  expressed  by 

yi^)  = 

where  e(z)  is  some  white  noise  sequence.  This  can  be  rewritten  as  in  Eq.  (5.2): 

MQ~^)yk  =  B(q~^)uk  +  A(q~^)ek 


However,  the  disturbance  A{q  ^)ek  is  no  longer  white.  To  “whiten”  the  system,  the  follow¬ 
ing  is  considered. 


Here,  the  input  and  output  are  filtered  through  1/A[q~^),  the  disturbance  is  once  again 
white,  and  a  LS  estimate  can  be  correctly  used.  The  problem  with  this  formulation  is  that 
A{q~^)  is  not  known.  To  overcome  this,  an  iterative  procedure  is  once  again  used.  The 
initial  estimate  is  simply  the  LS  estimate 


^^5  =  ^LS 


This  estimate  is  used  to  estimate  the  necessary  A{q~^)  and  the  input  and  output  are  then 
filtered  as  follows. 


yk  —  ®iyfc— 1  T  ®2yfe— 2  T  ■  ■  ■  T  dnyk—n  "t"  yk 
fifc  ~  1  T  ^2^k—2  T  ■  ■  ■  "h  O'n^k—n  T  '^k 

Next,  yk  and  Uk  replace  the  Y  and  U  in  the  H  matrix  to  form  the  H  matrix,  and  the  next 
estimate  is  calculated  as  follows: 

0gl5=(h''h)"'h^z 
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This  process  of  filtering  and  estimating  is  repeated  until  some  suitable  convergence  criteria 
is  met,  i.e.,  ||  0gl5  ~  ^GLS  II  sufficiently  small.  It  seems  that  this  method  also  con¬ 
verges  for  small  enough  noise  values,  provided  enough  measurements  are  available.  As  the 
noise  strength  is  increased,  more  and  more  measurements  are  needed  to  consistently  obtain 
an  estimate.  The  GLS  method  also  appears  to  suffer  from  the  problem  of  converging  to 
incorrect  estimates  for  low  enough  signal-to- noise  ratios. 


In  this  scenario,  it  is  assumed  that  there  is  noise  only  on  the  output,  and  that  the  input 
is  known.  Therefore,  according  to  the  GMV  paradigm,  the  equation  error  noise  vector  in 
Eq.  (5.6)  is  given  by 

h  =  Vk-  aiVk-i  -  a2Vk-2 
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Figure  5.1.  Discrete  Bode  Plot  of  Given  Plant 


To  obtain  the  equation  error  covariance  matrix,  one  calculates 


E  {0/;  — -j-}  —  \ 


(1  +  +  fll)  ) 

r  =  0, 

cr^(oi02  -  oi). 

r  =  1, 

(-02) , 

r  =  2, 

0, 

r  >2. 

(5.9) 


which  means  that  the  equation  error  covariance  R  is  the  pentadiagonal  matrix  R  =  a^'Ry, 
where 


R„  = 


a\  -ai  +  ai02  -a2 
—Cl  +  aia2  l-\-  a\-\-  a\  —a\  +  0102 
—02  — 1  +  «!  +  “2 


0 

-02 

—  Oi  +  O1O2 


0  0 
5.3  Algorithm  Comparison 


0 

0 

0 

1  +  +  ®2 


(5.10) 


In  this  section,  the  GLS  and  GMV  estimates  are  compared.  To  do  this,  Matlab’s 
[22]  randnO  function  is  initialized  with  a  seed  of  zero  and  used  to  generate  noise  with  a 
covariance  =  (10^  *  0.0001)^.  This  noise  is  then  added  to  the  true  output  obtained  from 
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the  system  given  in  Eq.  (5.8).  Five  levels  of  noise  {p  =  0,1, 2, 3, 4)  are  examined  in  this 
work,  but  only  three  are  presented  here.  The  low  noise  case  of  p  =  0  is  presented  first.  The 
output  for  this  case,  corrupted  by  representative  noise,  is  shown  in  Fig.  5.2. 

p=<) 


>- 


Time  (sec) 

Figure  5.2.  Plant  Output  with  Representative  Low  Noise  Level 


The  p  =  2  and  3  cases  are  not  presented  here  because  the  results  are  very  similar 
to  the  first  case.  Only  when  the  noise  level  increases  to  p  =  3  or  4  do  the  problems  and 
differences  become  apparent.  The  noise  levels  for  these  cases  are  shown  in  Fig.  5.3. 


Time  (sec)  Time  (sec) 

Figure  5.3.  Plant  Output  with  Representative  High  Noise  Levels 
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Figure  5.4  displays  the  denominator  estimates  for  a  100  run  Monte-Carlo  (MC) 
analysis  for  two  different  regression  lengths  -  40  and  70.  The  estimation  results  are  first 
normalized  by  dividing  each  estimate  by  the  true  estimate,  and  then  plotted.  Ellipses  are 
plotted  representing  each  estimation  method’s  actual  one  sigma  variation.  The  ellipses’  axes 
are  centered  at  the  average  estimate  for  each  method.  The  algorithm  predicted  estimation 
error  covariances  are  not  plotted,  because  they  are  close  to  the  actual  covariances  in  each 
case. 


p=0  m=40  p=0  m=70 


Figure  5.4.  Denominator  Estimates  for  40  (left)  and  70  (right)  Measurement  Linear  Re¬ 
gression 

As  can  be  seen,  the  GLS  estimates  are  quite  a  bit  better  than  the  GMV  estimates 
for  the  case  where  40  measurements  are  taken  (approximately  4  seconds  of  data).  However, 
this  does  not  teU  the  whole  story.  When  70  measurements  (approximately  7  sec  of  data)  are 
used,  both  estimates  are  better,  and  the  statistical  results  are  virtually  the  same.  To  better 
visualize  this  behavior,  Fig.  5.5  shows  the  GLS  and  GMV  estimation  error  and  covariance 
for  the  a\  parameter. 

As  can  be  seen,  the  GLS  estimate  is  initially  better  than  the  GMV  estimate.  How¬ 
ever,  within  about  60  measurements,  the  GMV  has  caught  up  with  the  GLS,  and  the  two 
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Number  of  Measurements  Number  of  Measurements 


Figure  5.5.  GMV  and  GLS  -  Incremental  a\  Estimation 

are  almost  identical  afterwards.  The  difference  is  even  smaller  when  the  initial  time  is  re¬ 
moved.  Figure  5.6  contains  the  incremental  estimation  of  the  ai  parameter  for  the  case  of 
no  initial  time. 


Number  of  Measurements  Number  of  Measurements 

Figure  5.6.  GMV  and  GLS  -  Incremental  ai  Estimation  without  Initial  Time 

From  this  information,  it  may  seem  that  the  GLS  is  always  at  least  as  good  as  the 
GMV,  but  the  cases  above  are  for  a  very  small  noise  level.  When  the  noise  level  increases, 
a  problem  appears  in  the  GLS  estimation.  The  GLS  method  requires,  in  general,  more 
measurements  initially  to  consistently  produce  a  valid  estimate.  Figure  5.7  contains  the 
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incremental  a\  estimation  error  and  covariance  for  the  case  of  p  =  3.  Once  again,  it  is 
apparent  that  the  GLS  estimate  is  better  than  the  GMV  up  until  about  60  measurements, 
but  what  is  not  shown  is  that,  while  the  GMV  algorithm  produces  an  estimate  for  m  <  10, 
the  GLS  algorithm  does  not  produce  an  estimate  for  aU  100  noise  realizations  until  m  =  19. 


GMV  -  al  Estimation  Error  (%) 
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Figure  5.7.  Incremental  oi  Estimation  for  p  =  3 


This  initial  measurement  requirement  is  more  visible  when  the  initial  time  is  removed 
as  well.  Figure  5.8  shows  these  results.  The  lack  of  information  in  the  GLS  plot  below  about 
m  =  25  is  due  to  the  fact  that  the  algorithm  does  not  produce  an  estimate  for  this  area. 
Similar  things  happen  for  the  case  of  p  =  4  as  well.  The  GLS  estimation  requires  even 
more  measurements  to  produce  an  estimate.  The  normalized  numerical  results  are  given  in 
Appendix  D  for  each  of  the  identification  algorithms  compared  in  this  work.  The  values  are 
for  p  =  3,4  with  400  measurements. 

Another  difficulty  with  both  methods  that  begins  to  play  a  major  role  at  this  noise 
level  is  the  problem  of  algorithm  convergence.  As  can  be  seen  in  Fig.  5.9,  the  algorithm 
predicted  covariance  is  quite  far  from  the  actual  covariance,  and  the  bias  is  rather  large 
below  60  measurements  for  both  algorithms.  This  is  due  to  the  fact  that  the  algorithm  did 
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GMV  -  al  Estimation  Error  (%)  GLS  -  al  Estimation  Error  (%) 


Number  of  Measurements  Number  of  Measurements 


Figure  5.8.  Incremental  aj  Estimation  for  p  =  3  without  Initial  Time 

not  converge  to  the  desired  fixed  point  when  initialized  with  the  LS  estimate.  Rather,  it 
converges  to  a  second  fixed  point,  as  is  discussed  in  the  next  section. 


GMV  -  al  Estimation  Error  (%)  GLS  -  al  Estimation  Error  (%) 


Number  of  Measurements  Number  of  Measurements 


Figure  5.9.  Incremental  a\  Estimation  for  p  =  4 

5.4  Convergence  Characteristics 

The  convergence  problem  in  the  time  domain  case  is  somewhat  worse  than  that 
in  the  frequency  domain.  In  the  frequency  domain,  there  was  one  small  area  around  the 
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origin  that,  if  the  algorithm  was  initialized  there,  caused  the  algorithm  to  converge  to  an 
incorrect  fixed  point.  The  solution  to  this  problem  was  just  to  not  initialize  in  that  area, 
but  rather  use  a  large  value  from  each  quadrant  and  choose  the  one  that  converged  the 
quickest.  However,  in  the  time  domain,  the  picture  appears  to  be  much  more  complicated. 
Figure  5.10  contains  the  estimation  results  for  m  —  40,  p  =  4.  As  can  be  seen,  both  the 
GMV  and  GLS  algorithms  suffer  from  this  problem,  and  the  erroneous  fixed  point  is  not 
necessarily  near  the  origin. 


p=4  m=40 


Figure  5.10.  Denominator  Estimates  Showing  Convergence  Problems 

There  are  at  least  three  different  time  domain  convergence  patterns.  One  is  the  case 
where  the  algorithm  converges  to  the  desired  estimate  no  matter  where  it  is  initialized.  The 
second  case  is  where  the  algorithm  converges  to  an  incorrect  fixed  point  when  initialized  at 
the  LS  estimate,  but  wiU  converge  to  the  desired  fixed  point  when  initialized  close  to  the 
true  parameter.  The  third  case  is  where  the  algorithm  wiU  not  converge  to  the  desired  fixed 
point  no  matter  where  it  is  initialized.  This  latter  case  is  truly  disconcerting,  and  it  is  not 
known  why  it  happens  or  how  to  correct  it.  The  only  things  that  can  be  said  about  it  are 
that  it  only  happens  for  a  small  number  of  noise  realizations  at  a  very  low  signal  to  noise 
ratio  where  there  are  not  yet  enough  measurements.  In  other  words,  there  is  insufficient 
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excitation.  The  phenomenon  disappears  once  “enough”  measurements  are  added  to  the 
linear  regression,  increasing  the  level  of  excitation. 

Two  examples  of  case  two  are  shown  in  Fig.  5.11.  They  are  the  6th  and  18th 
particular  noise  realizations.  The  different  areas  shown  in  the  plots  are  determined  using 
the  “shotgun”  approach  discussed  in  Section  3.6.3.  The  algorithm  is  initialized  at  different 
points,  and  the  thick  lines  seperate  the  areas  that  converge  to  different  points.  In  realization 
6,  the  small  area  around  the  true  parameter  is  the  area  that  yields  the  desired  parameter 
estimate  when  the  algorithm  is  initialized  there.  The  rest  of  the  space  converges  to  an 
incorrect  estimate.  In  realization  18,  the  upper  left  and  lower  right  quadrants  converge  to 
the  desired  estimate,  and  the  others  do  not. 


Realization  #6  Realization  #  1 8 


Figure  5.11.  Convergence  Pattern  for  Noise  Runs  6  and  18 

On  a  brighter  note,  looking  at  the  output  or  form  of  the  estimated  plants  may  shed 
some  light  on  the  multiple  fixed  point  quandry.  Figure  5.12  contains  the  true  output,  the 
measured  output  for  this  realization,  and  the  outputs  obtained  by  running  the  input  through 
the  two  different  estimated  plants.  As  can  be  seen,  the  correct  output  is  much  smoother 
initially  than  the  incorrect  one  in  both  cases. 


5-13 


Realization  #6 


Realization  #18 


Figure  5.12.  Outputs  for  Noise  Realizations  6  and  18 


The  high  measurement  noise  may  fool  the  identification  algorithm  into  thinking  that 
there  is  a  high-frequency  component  in  the  plant.  However,  upon  closer  inspection  of  noise 
realization  number  six,  it  can  be  seen  that  the  incorrectly  identified  plant  is  a  highly  unlikely 
candidate.  It  has  one  pole  on  the  positive  real  axis  and  one  on  the  negative  real  axis.  This 
is  not  a  possible  realization  for  a  physical  plant  with  a  zero-order  hold,  since  it  would  imply 
one  real  and  one  imaginary  pole.  Additionally,  in  this  case,  the  pole  on  the  negative  real 
axis  is  slightly  unstable. 


-0.94772  -I-  1.6921  -0.9477(z  -  1.7854) 

22  +  0.22902  -  0.7989  “  (2  4-  1.0156)(2  -  0.7867) 


(5.11) 


This  could  also  explain,  at  least  in  this  case,  why  adding  more  measurements  would  eliminate 
the  convergence  problems.  The  more  measurements  that  are  added,  the  less  the  unstable 
plant’s  output  would  match  the  measured  output. 

As  for  the  case  of  no  convergence  to  the  desired  parameter,  noise  realization  number 
20  is  a  good  example  of  this.  As  in  the  other  realizations,  there  are  two  fixed  points  present. 
Figure  5.13  shows  that  most  of  the  space  converges  to  an  unlikely  plant  similar  to  Eq.  (5.11), 
with  two  poles  on  the  real  axis  -  one  positive  and  one  negative. 


-1.85212 -f  3.2726  _  -1.8521(2-  1.7670) 

22  +  0.33892  -  0.5852  ~  {z  +  0.9530)(2  -  0.6141) 


(5.12) 
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Realization  #20 


Figure  5.13.  Convergence  Pattern  for  Noise  Run  20 


The  small  spot,  denoted  “other”,  converges  to  a  completely  different  plant,  but  one 
that  is  probably  as  unlikely  as  the  other.  This  identified  plant  is  given  by 


-5-95132  +  8.4410  -5.9513(2-1.4183) 

1^)206  -  ^2  _  0.18072  +  0.9685  ~  (2  -  0.0903  ±  iO.98) 


(5.13) 


As  can  be  seen,  the  imaginary  poles  are  now  both  plausible  for  a  zero-order  hold  plant,  but 
the  magnitude  of  the  poles  is  14.8  rad/sec.  If  a  plant  has  a  mode  with  a  frequency  this 
high,  then  0.1  sec  is  not  even  close  to  a  sufficiently  small  sampling  period  for  identification 
purposes.  The  output  plots  in  Fig.  5.14  bear  this  out  in  that  the  output  for  this  plant  looks 
very  similar  qualitatively  to  the  noisy  measured  output.  A  high  frequency  component  is 
readily  apparent. 

Therefore,  some  initial  knowledge/prior  information  is  necessary  if  one  is  trying  to 
identify  plants  with  a  very  small  number  of  measurements  for  the  noise  level  present.  If 
the  algorithm  fails,  it  does  not  seem  to  find  physically  believable  plants  for  the  assumed 
situation. 
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Realization  #20 


Figure  5.14.  Outputs  for  Noise  Realization  20 


5.5  Noise  on  Input 

An  additional  advantage  of  the  GMV  method  is  its  adaptability  to  other  linear 
regression  type  frameworks,  particularly  to  expanded  noise  models.  If  there  is  noise  on  the 
input  as  well  as  the  output,  then  the  GLS  derivation  is  completely  incorrect.  The  GMV 
algorithm  can,  however,  be  easily  adapted  to  handle  input  noise  simply  by  modifying  Vk  to 
include  the  new  noise.  Here, 

=  '^k  —  0,\Vk-l  —  a2Vk-2  —  biWk-1  —  b2'Wk-2 


Now,  the  needed  stochastic  modeling  entails  the  calculations 


E{VkVk-r}  =  < 


(1  -|-  af  -|-  af)  +  cr^,  (bj  +  bj)  , 

r  =  0, 

(aia2-  ai)  +  a^(bib2), 

r  =  1, 

r  =  2. 

.0, 

r  >  2. 

(5.14) 
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which  means  that  the  equation  error  covariance  matrix  R  is  now  the  pentadiagonal  matrix 
R  =  Ruj,  where  R„  is  given  in  Eq.  (5.10)  and  the  tridiagonal  R^y  is  as  follows: 

bl  +  bl  M2  0  •••  0 

6162  bj  +  bj  bib2  •••  0 

R^  = 

0  0  0  ■■■  bj  +  bl 

For  this  comparison,  the  strength  of  noise  on  the  input  is  the  same  as  that  on  the 

output.  For  the  case  of  very  small  input  noise,  the  GLS-GMV  comparison  is  similar  to  the 
case  of  no  input  noise.  The  GLS  is  better  for  small  m,  but  the  GMV  catches  up,  as  is  shown 
in  Fig.  5.15.  This  good  performance  by  GLS  is  probably  due  to  the  filtering  effect,  i.e.,  the 
effects  of  noise  can  be  reduced  by  filtering  the  input  and  output  through  low  pass  filters. 


(5.15) 


Figure  5.15.  Incremental  ai  Estimation  for  p  =  0  with  Input  Noise  and  No  Initial  Time 


However,  when  the  noise  level  increases,  the  input  noise  seems  to  have  a  more  pro¬ 
nounced  effect,  and  the  GLS  estimates  become  worse  than  the  GMV  estimates,  as  is  shown 
in  Figs.  5.16  and  5.17. 

Taken  to  the  extreme,  the  GLS  algorithm  appears  to  break  down.  Not  only  does  it 
fail  to  produce  an  estimate  for  small  m,  it  also  falls  for  large  m.  Even  when  it  does  produce 
an  estimate,  it  is  not  so  good  except,  strangely,  for  m  around  40  to  50.  The  GMV  estimate 


5-17 


Eiror  Mean  and  Covariance  Error  Mean  and  Covariance 


GMV  -  al  Estimation  Error  (%)  GLS  —  al  Estimation  Error  (%) 


Figure  5.16.  Incremental  a\  Estimation  for  p  =  3  with  Input  Noise 


GMV  -  al  Estimation  Error  (%) 


GLS  -  al  Estimation  Error  (%) 


—  MC  mean  error 
I  -  -  MC +/- sigma 

Algorithm  Computed  Sigma 


40  60  80 

Number  of  Measurements 


100 


Figure  5.17.  Incremental  Estimation  for  p  =  3  with  Input  Noise  and  No  Initial  Time 


is  having  trouble  with  convergence  issues  for  small  m,  more  so  with  the  included  initial 
time  case.  It  takes  between  100  and  200  measurements  for  the  initial  time  case  to  converge 
consistently  from  the  LS  estimate,  while  the  other  case  converges  consistently  for  60  to  70 
measurements. 


GMV  -  al  Estimation  Error  (%)  GLS  -  al  Estimation  Error  (%) 


Number  of  Measurements  Number  of  Measurements 

Figure  5.18.  Incremental  ai  Estimation  for  p  =  4  with  Input  Noise 


GMV  -  al  Estimation  Error  (%)  GLS  -  al  Estimation  Error  (%) 


Number  of  Measurements  Number  of  Measurements 

Figure  5.19.  Incremental  ai  Estimation  for  p  =  4  with  Input  Noise  and  No  Initial  Time 
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5. 6  The  Effects  of  Sampling  Rate 


In  this  section,  the  effects  of  sampling  rate  and  identification  time  interval  are  ex¬ 
amined.  To  do  this,  the  sampling  rate  for  the  system  is  increased  by  a  factor  of  4.  The  new 
discrete-time  transfer  function  is 


0.11922-0.1183 
~  22  -  1.97832-1-0.9792 


(5.16) 


The  incremental  identification  of  one  parameter  is  shown  in  Fig.  5.20.  The  left  side  contains 
the  identification  for  the  original  sampling  time  of  0.1  seconds,  and  the  right  hand  plot  is 
for  the  shortened  time  of  0.025  seconds.  As  can  be  seen,  the  number  of  measurements  in 
each  are  such  that  the  physical  time  interval  is  the  same.  It  is  obvious  from  the  plots  that 


GMV  -  al  Estimation  Error  (%) 


GMV  -  al  Estimation  Error  (%) 


Figure  5.20.  Incremental  ai  Estimation  for  0.1  (left)  and  0.025  (right)  Sampling  Times 


more  measurements  are  needed  to  obtain  the  same  quality  of  estimate  for  the  0.025  rate. 
Forty  measurements  at  0.1  produce  a  relatively  good  estimate,  while  the  estimate  at  forty 
measurements  for  0.025  is  quite  poor.  Therefore,  increasing  the  sampling  rate  will  worsen 
the  estimate  if  the  same  number  of  measurements  are  used. 

On  the  other  hand,  increasing  the  sampling  rate  while  retaining  the  same  physical 
identification  time  interval  produces  a  substantial  improvement  in  the  quality  of  the  es¬ 
timate.  Figure  5.21  contains  the  denominator  estimates  for  10  seconds  of  data  for  both 
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sampling  rates.  As  can  be  seen,  the  0.025  rate,  which  has  400  measurements,  is  approxi¬ 
mately  ten  times  better  than  the  0.1  rate,  which  has  only  100  measurements.  Conversely, 


p=3  m=100  p=3  in=400 


'j  yy-’t,  _ , _ i _  .  I  b _ , _ , _ I _ , _ , _ , — I 

0.998  0.999  1  1.001  1.002  0.9998  0.9999  1  1.0001  1.0002  1.0003 

al  al 

Figure  5.21.  Denominator  Estimates  for  10  Seconds  of  Data  -  0.1  left  and  0.025  right 


increasing  the  0.1  rate  estimate  to  400  measurements  only  results  in  a  factor  of  two  im¬ 
provement.  Therefore,  if  one  can  afford  the  computational  burden  of  the  increased  number 
of  measurements,  it  is  a  very  good  thing  to  have. 

5. 7  Conclusions 

In  this  chapter,  the  performance  of  the  Generalized  Minimum  Variance  and  Gen¬ 
eralized  Least  Squares  estimation  algorithms  is  compared.  Gaussian  measurement  noise  is 
assumed,  as  is  customary  in  statistical  filtering  and  system  identification  work.  The  novel 
minimum  variance  estimate  equations  are  derived  and  applied  to  a  nonlinear  estimation 
problem.  The  results  are  then  compared  to  the  established  Generalized  Least  Squares  esti¬ 
mate  for  a  second-order  system  that  is  representative  of  an  aircraft’s  pitch  dynamics,  which 
is  used  for  inner-loop  flight  control  system  design. 

While  the  GLS  method  requires  more  measurements  initially  to  obtain  an  estimate, 
it  appears  to  produce  a  better  estimate  at  this  point  than  the  GMV.  However,  the  two 
methods  appear  to  be  statistically  equivalent  as  more  measurements  are  added.  One  distinct 


5-21 


advantage  of  the  GMV  method  is  its  ability  to  produce  valid  estimates  at  small  measurement 
samples,  and  in  general,  under  conditions  of  poorer  excitation. 

The  GMV  and  GLS  estimates  outperformed  the  LS  estimate  in  all  cases.  The  LS 
estimate  did,  however,  provide  a  useful,  albeit  sometimes  dangerous,  starting  point  for  iter¬ 
ating  the  parameter  estimate.  At  small  noise  levels,  it  does  not  matter  where  the  algorithms 
are  initialized,  but  as  noise  levels  increase,  there  is  an  additional  point  of  convergence  that 
can  trap  the  estimate.  In  high  noise  cases  (low  signal-to- noise  ratio),  it  becomes  necessary 
to  examine  the  final  result  to  see  if  it  is  a  valid  estimate.  There  seem  to  be  no  definitive 
differences  between  the  convergence  rates  for  incorrect  verses  correct  parameter  estimates 
in  the  time-domain.  Obviously,  there  must  be  at  least  a  small  amount  of  prior  information 
about  the  plant  for  a  valid  estimate  to  be  obtained. 

As  before,  the  least  squares  estimate  is  not  effective,  even  in  cases  of  small  mea¬ 
surement  noise.  The  GLS  and  GMV  provide  much  more  accurate  estimates,  and  the  GMV 
seems  to  provide  estimates  with  much  fewer  measurements.  Additionally,  the  GMV  is  much 
more  adaptable  to  various  identification  requirements. 
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VI.  Identification  of  the  Physical  Parameters  of  a  Continuous-Time  Dynamical 

System 

In  this  chapter,  the  identification  of  physical  parameters  using  the  sampled  input 
and  output  of  continuous-time  plants  is  discussed.  The  experiment  is  based  on  Case  3  of  a 
system  identification  competition  out  of  Italy  [3],  and  is  related  to  the  work  in  [5].  These 
consider  the  identification  of  thermal  building  parameters,  which  are  typically  modeled  as 
a  second-order  thermal  circuit  with  two  inputs  and  one  output. 

Currently,  an  important  problem  of  classical  system  identification  is  the  methods  by 
which  algorithms  are  validated.  Too  often,  the  researchers  start  with  a  set  of  input /output 
data  from  an  unknown  plant.  They  use  the  data  to  estimate  the  parameters  of  the  plant,  and 
then  use  the  input  data  and  the  identified  plant  to  simulate  the  output.  If  the  original  output 
and  the  simulated  output  are  “close” ,  then  the  algorithm  is  considered  to  be  successful.  The 
flaw  in  this  validation  “method”  becomes  apparent  under  conditions  of  less  than  optimal 
excitation.  With  poor  excitation,  it  is  possible  for  drastically  different  plants  to  produce 
the  same  output  from  the  same  input.  This  is  a  typical  manifestation  of  an  inverse  problem, 
which  system  identification  is. 

The  proposed  proper  method  of  algorithm  validation  is  to  start  with  a  known  plant, 
produce  an  input /output  data  set  using  this  plant,  add  measurement  noise  commensurate 
with  the  known  prevailing  signal-to-noise  ratio  in  the  actual  experiment,  identify  the  pa¬ 
rameters,  and  compare  these  parameters  to  the  original  known  plant.  If  the  parameter 
estimates  are  good  in  this  “synthetic”  experiment,  then  the  algorithm  can  be  applied  to  the 
unknown  plant  as  well. 

Section  6.1  is  a  summary  of  the  problem.  Section  6.2  explains  the  method  used  for 
identification  of  the  discrete  time  parameters.  Section  6.3  explains  the  procedure  of  the 
experiment.  Section  6.4  explains  how  to  use  the  initial  knowledge  of  the  system  to  improve 
the  identification,  and  Section  6.5  contains  the  methods  used  for  obtaining  the  physical 
parameter  estimates. 
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6.1  Physical  System 


The  physical  system  for  this  experiment  is  represented  by  the  second-order  linear 
thermal  network  shown  in  Fig.  6.1.  Ti  and  Tg  are  the  internal  and  external  temperatures, 


Figure  6.1.  Layout  of  the  Physical  System 


respectively,  Q  is  the  heat  flow,  iZi,  R2,  and  R3  are  thermal  resistance  elements,  and  Gi 
and  G2  are  thermal  capacitance  elements.  In  the  competition,  each  data  set  consisted  of  30 
days  of  hourly  observations  of  temperature  and  heat  flow.  A  relatively  high  level  of  noise 
was  added  to  Tg,  a  low  noise  level  was  added  to  Q,  and  no  noise  was  added  to  Ti.  The 
noise  on  Tg  results  in  a  variation  of  the  usual  system  identification  problem,  where  the  input 
signals  are  assumed  known  (noiseless). 

The  two  states  in  this  system  are  Ti  and  T2,  the  two  inputs  are  Tg  and  Ti,  and  the 
output  is  Q.  Therefore,  letting  Hi  =  1/Ri,  i  =  1,2,3,  the  state  space  equations  are 

Gi^  =  HiiTe  -  Ti)  -1  H2{T2  -  Ti) 

G2^  =  H2{Tr  -  T2)  +  H^{Ti  -  T2) 

Q  =  H^{Ti  -  T2) 


or  in  the  more  familiar  form  x  =  Ax  -f  Bu  and  y  =  Cx  -1-  Du, 


1 

r  Hj+Hj 
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(6.1a) 


(6.1b) 
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This  continuous-time  system  is  a  Two-Input  Single-Output  (TISO)  system  that  is 
not  strictly  proper.  The  state  equations  given  above  yield  the  following  continuous-time 
transfer  function: 


Q{s) 


Bqs^  +  Bis  -f  1 
AqS^  -f  Axs  -1-  ^2  * 


(6.2) 


where 


Bq  =  RiR^GiG^  (6.3a) 

=  iZiGi  +  R2G2  +  R1G2  (6.3b) 

Aq  =  R\R2R^G\G2  (6.3c) 

Ai  =  R1R2G1  -|-  RiR^Gi  -|-  R1R3G2  -f  R2R3G2  (6.3d) 

A2  —  R\  A  R2  T  R3  (6.3e) 


The  inputs  are  Ti  and  Tg,  the  output  is  Q,  and  the  five  unknown  physical  parameters  to  be 
identified  are  ili,  R2,  R3,  Gi,  and  G2. 


6.2  Time  Domain  Estimation  Method 


When  a  Tustin  transformation  is  performed  on  a  proper  or  strictly  proper  second- 
order  continuous  system,  such  as  that  in  Eq.  (6.2),  the  resulting  discrete  time  system  is  not 
strictly  proper.  In  other  words,  the  following  discrete  system  is  the  result. 


9(^) 


bipz'^  A  h\iz  -|-  612 
—  a\Z  —  02 


ti{z)  A 


^20^^  A  b2iz  A  ^22 
z^  —  a\z  —  02 


te{z) 


(6.4) 


Here,  the  discrete-time  transfer  function’s  8  parameters  oi,  02,  610,  &11,  bi2,  62O)  &21?  and 
622  are  unknown.  Although  the  original  continuous-time  system  had  5  parameters,  the 
Tustin  transformation  produces  8  in  the  discrete-time  system.  This  transfer  function  can 
be  rewritten  in  recursive  form  as 


Qk  —  +  0'2qk-2  +  bwtik  A  b\\tik-l  -|-  bi2tik-2  +  &2ofeA:  +  &2liei:-l  +  &224fe-2 


(6.5) 
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If  one  can  obtain  the  true  output  and  the  true  inputs  tik  and  tek,  then  only  eight 
measurements  are  needed  to  solve  for  the  unknown  parameters.  However,  the  true  values 
cannot  be  directly  measured.  Rather,  the  noise  corrupted  Qk,  Tik,  and  Tek  are  measured, 
where  the  noise  is  assumed  temporally  uncorrelated  (white)  and  Gaussian: 

Qk  =  qk  +  Vk,  V  =  j\r(o,  al) 

Tik  —  ^ik  "t"  '^ik, 

Tek  —  ^ek  ‘^ekt  '^e  — 

which,  when  substituted  into  Eq.  (6.5),  yields 

Qk  ^k  —  ®l(Qfc— 1  l)  ®2(Qa;— 2  2)  "b  '^ik^  "b  Wik—l') 

+  bi2{Tik-'2  —  Wik-2)  +  b2o{Tek  —  Wek)  +  62l(T’efc-l  ~  '^ek-l)  +  b22{Tek-2  ~  Wek-2) 


Qk  —  (l‘lQk-l+0,2Qk-2+bioTik  +  biiTik-i  +  bi2Tik-2  +  b2oTek  +  b2lTek-l+b22Tek-2  +  ^k  (6-6) 

where 

H  =  ‘^k-  O.lVk-1  ~  a2Vk-2  -  b\QWik  -  biiWik-i  -  bi2'Wik-2  -  b20Wek  ~  b2lWek~l  “  b22Wek-2 

This  expression  is  now  in  a  form  that  can  be  set  up  in  a  statistical  linear  regression 
equation  given  by  Eq.  (2.7),  where 

Qk 

Qk+l 

Z  = 

Qk+N 


Qk~l 

Qk-2 

Tik 

Tik-i 

Tik-2 

Tek 

Tek-1 

Tek-2 

H  = 

Qk 

Qk-1 

Tik+i 

Tik 

Tik-\ 

Tek+1 

bs”  ••• 

Tek—l 

Qk-l+N 

Qk-2+N 

Tik+N 

Tik-l+N 

Tik-2+N 

Tek+N 

Tek—l+N 

Tek-2+N 

6-4 


-iT 

^  =  ai  a2  bio  bii  612  620  ^>21  &22 

h 

4+1 

€  = 

i’fc+iV 

The  R  matrix  is  the  expected  value  of  ee^,  and  if  the  noise  is  assumed  white,  then 
the  R  matrix  is  an  iV  x  iV  Toeplitz  matrix  with  a  non-zero  diagonal  and  2  non-zero  off- 
diagonal  terms  above  and  below  the  diagonal.  If  v  and  We  are  assumed  uncorrelated  and 
the  input  measurement  Ti  is  noiseless  =  0),  then  the  R  matrix  is  given  by 

1  +  al  +  a2  —ai  +  aia2  — a2  •••  0 

—  til  T  ®1®2  1  d"  d”  ®2  — d"  tliti2  •  ■  •  0 

-02  -Old- 0102  1 -b  of -f  O2  •••  0 

0  0  0  •  •  •  1  d-  d-  ttf 

f>20  d-  ^>21  d-  ^>22  ^20^21  d"  ^11^21  &20&22  ‘  •  0 

^20^21  d-  ^*11^21  ^20  d-  &21  d-  ^>22  d"  1>11&21  ‘  ‘  0 

d"^u)e  f>20^22  ^20^21  d"  &11&21  ^>20  d"  ^»21  d"  ^>22  '  ’  ’  ^ 

0  0  0  •  •  •  620  d"  ^*21  d"  ^*22 

6.3  Problem  Specifics 

The  inputs  used  in  this  problem  are  created  by  taking  the  first  five  or  six  peak 
frequencies  from  a  Fourier  analysis  of  the  competition  inputs.  Therefore,  the  inputs  used 
in  this  experiment  are  just  a  sum  of  cosines  and  are  given  by 

6 

Tqi  —  ^  Moik  cos  {uJakl  T  ) ?  CX  =  i,  6 
k=l 

where  M,  w,  and  (j)  are  given  in  Table  6.1. 
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Table  6.1.  Values  for  the  Ti  and  Te  Inputs 


Q 

Mik 

>^ik 

^ik 

Alek 

4^ek 

3 

25.4132 

0 

0 

11.0476 

0 

0 

IQI 

-0.9391 

0.0245 

-2.0109 

0.0123 

-0.9026 

ill 

-0.7854 

0.0368 

-0.8515 

0.0245 

-2.3968 

-0.8098 

0.2454 

-3.0426 

-2.7772 

0.2577 

-3.0037 

Ql 

-2.3719 

0.2577 

2.9906 

-1.6586 

0.2700 

0.6751 

Ql 

-1.2646 

0.2700 

0.3943 

0 

0 

0 

These  inputs  are  given  to  the  system  described  in  Eq.  (6.2),  where  the  physical 
parameters  are  those  of  the  competition:  R\  —  1,  R2  =  0.1,  R^  =  10,  G\  =  100,  and 
G2  =  50.  Next,  1000  periods  of  data  are  taken,  but  only  the  final  402  (m  =  400)  points  are 
used.  This  is  done  in  an  attempt  to  simulate  a  building  that  is  in  steady  state. 

Finally,  a  low  noise  level  is  added  to  the  output,  and  a  higher  noise  level  is  added 
to  the  external  temperature  input.  The  signal  to  noise  ratios  are  based  on  the  competition, 
and  are  chosen  to  be  130  dB  and  50  dB,  respectively. 

6.4  A  Priori  Information 

If  one  attempts  to  “brute  force”  identify  aU  eight  of  the  parameters  in  Eq.  (6.4), 
then  the  Te  numerator  results  are  very  poor.  This  is  probably  because  the  discrete  system 
is  over-determined.  There  are,  in  actuality,  only  five  parameters  in  the  system.  Keeping  this 
in  mind,  one  can  modify  the  system  as  follows  so  only  five  parameters  need  be  identified. 

6.4- i  Modification  of  the  Known  Continuous-Time  System.  Rather  than  using 
brute  force,  it  is  useful  to  use  the  a  prior  knowledge  that  one  has.  To  this  end,  rewrite 
the  transfer  functions  in  Eq.  (6.2)  as 

(BqsH  Bis)  r,(s)  Ti{s)~Te{s) 

^  Aos2  -f  Ais  +  A2  Aos2  -f  Ais  +  A2  ^  ^ 
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When  a  Tustin  transformation  is  performed  on  this  system,  the  following  system  is  the 
result. 

=  "72— 77— +  ;2  „  ,  „ —  -  ^-(^1  6.8) 

—  a,\z  —  02  z‘‘  —  a\z  —  02 

It  would  appear  that  there  are  still  eight  parameters  to  identify,  but  when  a  Tustin  trans¬ 
formation  is  performed  on  a  second-order  continuous-time  transfer  function  with  a  constant 
numerator,  as  in  the  second  term  of  Eq.  (6.8),  621  =  2&20  and  622  =  ^20-  Also,  when  a  Tustin 
transformation  is  performed  on  a  second-order  numerator  with  the  form  Ks{s  -f  z),  as  in 
the  first  term  of  Eq.  (6.8),  611  +  612  +  613  =  0.  Therefore,  Eq.  (6.8)  is  really 

/  s.  i—bn~bi2)z^  +  bnz  +  bi2  ,  .  b2oz^  +  2b2oz  +  b2o  .  ,  .  ^  rw 

^ - o  '  . - +  -  .2  -  -  ; - (6-9) 


Z^  —  OiZ  —  02 


Z^  —  0\Z  —  O2 


which  changes  Eq.  (6.6)  to 


where 


Qk  —  O.lQk-1  +  0,2Qk-2  +  bi\{Tik-\  —  Tik)  +  bi2(Tik-2  —  Tik) 

+b20  {{Tik  -  Tek)  +  2{Tik-\  —  Tek-l)  +  {Tik-2  -  Tek-2))  +  ^k  (6.10) 


^k  —  '>^k  ~  ~  02^^fc-2  ~  {~btl  —  bi2  -f  b2o)Wik  ~  {bn  +  2b2o)Wik-l 

—  {bl2  +  b2o)Wik-2  +  b20Wek  +  2b2oWek-l  +  ^20'l<^efc-2 


The  number  of  discrete  parameters  to  be  estimated  has  been  decreased  to  five,  and  is  now 
equal  to  the  number  of  unknown  physical  parameters  in  the  original  system.  This  changes 
the  H  and  R  matrices  to  {a^j^i  =  0  and  T^k  =  life  -  Tek) 

r  Qk-l  Qk-2  {Tik-1  -  Tik)  {Tik-2  -  Tik)  {TAk  +  ^TAk-1  +  TAk-2)  1 


R  =  cr2 


1  -f-  -f-  (I2  T  ®ifl2  —02 

— Oi  -(-  0\02  1  -f  uj  O2  — fll  ~\~  O1O2 

-O2  — +  O1O2  o\-\-  02 


1  +  Oj  -t- 
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6  4  1  •••  0 

4  6  4  •••  0 

1  4  6  •••  0 
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6.4-2  Feasibility  of  the  Tustin  Transformation.  Of  concern  in  this  work  is  the  va¬ 
lidity  of  the  Tustin  transformation.  The  Tustin  transformation  assumes  ramp  type  inputs, 
and  substitutes  a  bilinear  approximation  for  exp  (At)  in  the  conversion  process  (see  Ap¬ 
pendix  A).  The  inputs  to  the  discrete  Tustin  system  would  be  the  same  as  the  clean  sampled 
inputs  to  the  continuous  system.  Therefore,  the  difference  would  be  in  the  sampled  out¬ 
puts.  To  try  to  determine  this  difference,  the  discrete  system  in  Eq.  (6.9)  is  simulated  with  a 
sampled  version  of  the  continuous  input  that  is  used  with  the  continuous-time  system.  The 
results  indicate  that  there  is  a  discrepancy  between  the  discrete  and  continuous  outputs  as 
shown  in  Fig.  6.2.  As  is  seen,  there  is  a  relatively  large  initial  difference  between  the  two 


Figure  6.2.  Comparison  of  Continuous  and  Discrete  Outputs 

system  outputs  if  the  entire  data  set  is  used.  This  agrees  with  results  discussed  in  Ref.  [13]. 
This  difference  is  probably  due  to  the  initial  transient  in  the  continuous  system,  which  is 
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too  fast  for  the  Tustin  transformation  to  model.  The  large  initial  difference  introduces  a 
relatively  large  error  in  the  estimate,  and  it  takes  a  very  large  estimation  window  to  wash 
it  out. 

However,  if  the  system  measurements  are  taken  after  the  system  is  in  steady  state, 
the  initial  differences  are  smaller.  If  the  corresponding  discrete  system  is  simulated  using 
only  the  information  in  the  final  402  measurements,  then  the  results  in  Fig.  6.3  are  obtained. 


Figure  6.3.  Difference  with  Initial  Time  Removed 

To  examine  the  possible  effects  of  these  continuous/discrete  differences,  the  identifi¬ 
cation  experiment  is  run  with  a  very  small  measurement  noise  on  both  the  Tg  input  and  the 
output,  for  both  the  continuous  and  corresponding  discrete  systems.  The  identification  re¬ 
sults  for  the  continuous  system  are  shown  in  Fig.  6.4,  and  the  discrete  results  are  in  Fig.  6.5. 

As  can  be  seen,  the  induced  biases  are  much  less  than  1%  ,  so  the  Tustin  transformation 
is  considered  satisfactory  for  this  problem.  If  the  sampling  rate  were  lower,  or  if  the  initial 
time  information  was  needed,  then  perhaps  a  different  continuous-discrete  transformation 
would  be  needed.  Other  noise  models  are  investigated  in  this  work  as  well  in  an  attempt  to 
reduce  the  error. 
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It  can  be  seen  in  Fig.  6.3  that  the  difference  between  the  continuous  and  discrete 
outputs  consists  mainly  of  a  high  frequency  sinusoid  and  also  possibly  of  a  ramp  in  the 
beginning.  There  are  at  least  three  possible  ways  to  treat  this.  One  could  model  the 
disturbance  as  a  sinusoidal  error,  a  combination  sinusoid/ramp,  or  simply  as  white  noise. 
Each  of  the  three  were  added  to  the  system  model,  but  the  extra  parameters  added  by  the 
additional  modeling  seem  to  negate  any  gains,  and  little,  if  any,  improvement  is  achieved. 
Therefore,  in  order  to  keep  the  number  of  estimated  parameters  to  a  minimum,  a  simple 
white  noise  model  is  used. 

6. 5  Estimation  of  the  Physical  System  Parameters 

This  section  discusses  the  procedure  used  to  obtain  the  required  estimates.  Once 
the  discrete  system  parameters  and  error  covariances  are  estimated,  they  are  converted  into 
the  estimates  and  covariances  for  Gi,  G2,  and  R^.  To  accomplish  this,  the  discrete 

parameters  are  first  converted  to  their  continuous-time  equivalents  via  an  inverse  Tustin 
transformation,  and  then  a  set  of  nonlinear  equations  is  solved  for  the  physical  parameters. 
The  covariance  estimates  are  obtained  via  perturbation  methods. 

6.5.1  Identification  of  the  Parameters  of  the  Discretized  System.  As  discussed  in 
Section  6.4.2,  only  the  final  402  points  are  used  from  the  data  sets  to  simulate  a  building 
in  steady-state.  Next,  a  100  run  Monte  Carlo  analysis  is  performed.  For  each,  the  least 
squares  estimate  is  calculated  as  an  initial  guess,  and  the  Generalized  Minimum  Variance 
estimate  is  iterated  fifty  times.  Figure  6.6  shows  the  average  estimate,  each  of  the  100 
estimates,  and  the  calculated  covariance  for  both  of  the  denominator  parameters  and  two 
of  the  numerator  parameters.  Table  6.2  contains  the  normalized  numerical  results  for  the 
bias  (e),  actual  error  covariance  (ce),  and  algorithm  predicted  covariance  (ap)  for  all  the 
parameters. 
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6.5.2  Conversion  to  Continuous-Time.  The  first  order  of  business  here  is  to  use 


an  inverse  Tustin  transformation  to  convert  the  discrete  parameters  of  Eq.  (6.9)  into  a  set 
of  continuous-time  parameters  as  in  Eq.  (6.2). 

The  calculation  for  the  error  covariance  estimate  of  the  continuous-time  param¬ 
eters  is  done  using  an  approximation  of  the  Jacobian  matrix  J.  To  do  this,  let  9c  = 
[Ao  Ai  A2  Bo  Bi]^  and  Oj,  =  [ui  bn  bn  &2o]^-  Then  J  is  numerically  calculated  such 
that  AOc  =  J  AOd,  where  =  (^ed[i]-  This  approximation  is  exact  for  linear  systems, 

and  it  should  be  close  for  slightly  nonlinear  systems.  It  is  assumed  that  the  systems  here  are 
slightly  nonlinear  for  the  purpose  of  Jacobian  matrix  computation.  After  these  calculations 
are  performed,  the  covariance  of  the  continuous  parameters  is  given  by  Pc  =  JP^J^.  The 
normalized  continuous  parameter  estimates  and  covariances  are  shown  in  Fig.  6.7,  and  the 
normalized  numerical  results  are  in  Table  6.3. 

m=400  m=400 


Figure  6.7.  Continuous  Parameter  Estimates  and  Error  Covariances 


6.5.3  Conversion  to  Physical  System  Parameters.  After  the  estimate  and  covari¬ 
ance  for  6c  is  obtained,  the  set  of  nonlinear  equations  given  in  Eq.  (6.3)  are  solved  for  the 
physical  system  parameters. 


R3  = 


Aq 

Bo 
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Table  6.3.  Normalized  Numerical  Results  for  Continuous  Parameters 


e 

(^e 

(Tp 

^0 

-0.3982 

2.1402 

1.9909 

ESI 

0.0054 

0.6181 

0.6456 

0.0044 

0.0571 

0.0647 

Bo 

-0.3756 

2.1501 

2.0015 

Bl 

-0.0022 

0.6326 

0.6599 

C2  = 
R2  = 


AiBq  —  AqB\ 

Bo 

C2B1  -  CiiA2  -  Rs) 
Ri  =  A2  —  R2  —  R3 
Bo 


Ci  = 


R1R2C2 


Once  the  parameter  estimates  are  obtained,  the  error  covariances  are  estimated 
similar  to  the  continuous  parameter  error  covariance  estimates,  i.e.,  A9  =  J  A^c,  where 
A0c[i]  =  <yec[i]-  The  results  of  these  final  calculations  are  shown  in  Fig.  6.8,  and  the 
numerical  results  are  in  Table  6.4. 


m=400  m=400 


Figure  6.8.  System  Parameter  Estimates  and  Error  Covariances 


6-14 


Table  6.4.  Normalized  Numerical  Results  for  Physical  Parameters 


e 

e  =  6~e{%) 

o-e 

Op 

Ri 

0.0479 

0.6549 

0.9461 

R2 

2.2519 

5.5391 

5.6649 

Rs 

-0.0224 

0.0262 

0.0262 

2.1800 

4.3873 

4.1068 

(J2 

-4.2527 

6.7806 

5.2851 

6. 6  Conclusions 


This  chapter  investigates  a  method  of  identifying  physical  parameters  using  sampled 
measurements  of  continuous-time  inputs  and  outputs.  The  flexibility  of  the  GMV  algorithm 
is  apparent  in  this  problem,  adapting  easily  to  measurement  noise  on  one  input  and  on  the 
output.  Also,  the  system  has  been  rearranged  such  that  only  live  parameters  need  be  iden- 
tifled.  The  well  known  Tustin  transformation  is  used  to  convert  the  estimated  discrete-time 
parameters  to  continuous-time  parameters.  The  error  covariances  are  estimated  numeri¬ 
cally,  and  the  results  are  close  to  the  actual  covariances.  The  physical  parameters  and  error 
covariances  are  calculated  in  a  similar  manner. 

An  important  point  of  this  chapter  is  the  number  of  parameters  to  be  identified. 
If  one  tries  to  identify  an  over-determined  system,  the  results  wiU  suffer,  if  they  can  be 
obtained  at  all. 
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VII.  Conclusions 


This  research  identifies  and  focuses  on  the  solution  of  the  problems  that  measure¬ 
ment  noise  introduces  to  the  identification  of  dynamical  systems.  Measurement  noise  is  a 
fact  of  life  in  experimental  and  applied  science,  and  should  be  specifically  addressed.  The 
Generalized  Minimum  Variance  (GMV)  identification  procedure  developed  in  this  work  is 
one  possible  solution.  An  advantage  of  this  method  is  the  fact  that  it  is  based  on  linear 
regression  and  linear  prediction.  This  makes  it  easy  to  understand  and  adapt  to  other  es¬ 
timation/system  identification  problems,  where  it  produces  equally  good  results.  Further 
applications  of  the  GMV  algorithm  can  be  found  in  Ref.  [15]. 

This  work  introduces  the  concept  that  the  parameter  estimate  provided  by  the  GMV 
algorithm  is  a  fixed  point  of  a  nonlinear  mapping  derived  from  the  theory  of  Minimum 
Variance  identification.  Existence  of  a  fixed  point  is  proven,  and  the  particular  algorithm 
used  corresponds  to  the  iteration  suggested  by  a  contraction  mapping.  This  eliminates  the 
need  for  global  optimization,  and  the  problems  of  multiple  local  minima  are  greatly  reduced. 
In  fact,  only  one  other  possible  fixed  point  solution  was  observed  in  the  course  of  this  work 
at  low  signal-to-noise  ratios.  However,  this  benign  situation  may  not  persist  for  higher  order 
systems.  Moreover,  when  the  alternative  solution  is  found,  it  is  not  a  plausible  one  for  the 
given  situation,  and  possible  means  of  eliminating  it  are  suggested. 

In  Chapter  III,  the  GMV  algorithm  is  applied  to  the  identification  of  continuous¬ 
time  transfer  functions  using  frequency  domain  measurements.  Three  different  noise  models 
are  examined,  and  the  results  are  compared  to  those  from  the  classical  Least  Squares, 
Generalized  Least  Squares  (GLS),  and  Instrumental  Variable  methods.  Only  the  GLS 
algorithm  produces  similar  quality  parameter  estimates.  Moreover,  the  way  in  which  the 
GLS  method  is  cost  minimization-based  perhaps  clouds  the  treatment  of  measurement  noise, 
and  the  results  suffer.  When  the  noise  is  properly  treated  in  the  GMV  method,  the  results 
are  better.  Convergence  rates  for  the  two  algorithms  are  similar,  with  a  small  number  of 
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iterations  needed  at  small  noise  levels.  However,  large  numbers  of  iterations  are  needed  as 
the  noise  level  becomes  high. 

Prior  knowledge  of  the  system  order  is  needed  to  use  the  GMV  method,  and  some 
knowledge  of  the  system  bandwidth  is  needed  to  obtain  the  parameter  estimates  accurately. 
As  is  seen  in  Chapter  IV,  if  ones  takes  measurements  outside  the  bandwidth  of  interest, 
totally  inaccurate  estimates  are  obtained. 

In  Chapter  V,  the  GMV  algorithm  is  applied  to  the  identification  of  a  discrete-time 
plant.  This  scenario  is  applicable  to  control  system  design,  since  the  input  to  a  computer 
controlled  plant  is  typically  zero-order  hold  in  nature,  and  thus  one  obtains  a  direct  and 
correct  transformation  from  the  continuous-time  model  to  the  discrete-time  model.  The 
versatility  of  the  GMV  identification  method  is  apparent  in  this  situation  as  it  easily  adapts 
to  accommodate  the  addition  of  measurement  noise  on  the  input. 

Also  readily  visible  in  this  chapter  is  the  importance  of  correct  model  validation.  In 
general,  simply  because  the  estimated  and  measured  outputs  are  close  does  not  mean  that 
the  system  has  been  correctly  identified.  In  the  GMV  estimation  method,  an  iteration  is 
performed  to  find  a  fixed  point.  As  the  noise  increases,  another  fixed  point  appears.  When 
the  output  of  this  second  estimated  plant  is  compared  to  the  correctly  estimated  plant,  little 
quantitative  difference  is  apparent.  Thus,  the  output  of  two  distinctly  different  plants  can 
have  similar  outputs  for  a  given  input.  Fortunately,  at  least  in  this  scenario,  the  incorrectly 
identified  plant  is  not  a  plausible  solution. 

Finally,  Chapter  VI  applies  the  GMV  algorithm  to  the  identification  of  a  continuous¬ 
time  dynamical  system,  and  to  the  estimation  of  its  physical  parameters.  Here,  noise  is 
present  on  one  of  the  inputs  and  the  output  of  the  second-order,  two  input-one  output 
plant.  The  Tustin  transformation  is  determined  to  be  a  suitable  one  in  this  case,  because 
the  deterministic  errors  introduced  are  less  than  one  percent.  Note,  one  should  not  simply 
apply  any  identification  method  blindly.  If  only  actual  measurements  are  available,  with  no 
idea  of  the  actual  parameters,  one  should  first  perform  a  synthetic  experiment  and  gauge 
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the  results.  If  problems  arise  in  such  an  experiment,  they  can  then  be  addressed  before  one 
attempts  to  identify  the  actual  unknown  parameters  using  real  data. 

An  additional  important  aspect  of  system  identification  is  also  addressed  in  this 
chapter:  the  concept  of  keeping  the  number  of  identified  parameters  to  a  minimum.  This 
is  where  prior  knowledge  of  the  system  is  useful.  In  the  problem  addressed  in  Chapter  V, 
there  were  actually  only  five  parameters  describing  the  system.  If  the  problem  is  addressed 
in  the  time  domain  with  no  prior  knowledge,  then  the  number  of  needed  parameters  is  eight. 
Using  knowledge  of  the  system  however,  one  can  reduce  the  number  of  needed  parameters 
in  the  time  domain  to  five  -  the  minimum  number.  Failing  to  do  so  severely  hampers  the 
identification  process. 

In  conclusion,  the  attributes  of  proper  (stochastic)  modeling  cannot  be  overempha¬ 
sized.  The  performance  of  the  GMV  system  identification  algorithm  is  shown  to  be  equal  to 
or  superior  than  the  best  available  linear  regression  based  system  identification  methods.  It 
produces  good  parameter  estimates  for  small  numbers  of  measurements,  along  with  an  ac¬ 
curate  prediction  of  the  estimation  error,  when  convergence  is  not  a  problem.  Convergence 
of  the  GMV  algorithm  occurs  for  relatively  low  signal-to-noise  ratios.  Hence,  the  GMV 
algorithm  is  completely  autonomous  for  reasonable  noise  levels.  Also,  the  structure  of  the 
R  matrix  may  allow  the  design  of  recursive  algorithms  if  so  desired.  Finally,  being  based 
on  proper  stochastic  modeling,  the  GMV  system  identification  method  developed  in  this 
dissertation  is  shown  to  be  readily  adaptable  and  equally  applicable  to  a  whole  spectrum 
of  identification  and  estimation  problems. 
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Appendix  A.  Methods  for  Discretizing  Continuous  Systems 


While  physical  systems  are  almost  always  continuous-time  systems,  they  are  typically 
converted  into  discrete-time  models  for  many  purposes,  including  digital  control  and  system 
identification.  There  are  many  variations  in  discretization  methods,  and  most  are  based  on 
differing  assumptions  about  the  form  of  the  input. 

A.l  Continuous- Time  Model 

All  linear  (or  linearized)  physical  systems  can  be  represented  by  the  familiar  form 

x(t)  =  Ax(t)  +  Bu(t)  (A. la) 

y(t)  =  Cx(t)  -1-  Du(t)  (A.lb) 

where  x  is  a  vector  of  the  states  in  the  system,  u  is  a  vector  of  the  inputs  to  the  system, 
and  y  is  a  vector  containing  the  outputs  of  the  system.  The  goal  of  discretizing  the  system 
is  to  obtain  a  model  of  the  form. 

Xfc  =  AdXfc_i  +  BdUk-i  (A.2a) 

yfc  =  CdXfc  -f-  DrfUfc  (A.2b) 

where  X/t  =  'x.{kT),  Xfe_i  =  x((fc  -  l)r),  and  likewise  for  Ufc  and  y^.  Comparing  Eq.  (A. 2b) 
to  Eq.  (A.lb),  it  can  be  seen  that 

Cd  =  C  and  Dd  =  D  (A.3) 

To  solve  for  the  time  history  of  the  output,  given  the  time  history  of  the  input  u(t) 
and  the  system  states  at  the  initial  time  xq  =  x(to)5  the  following  is  used  [7]: 

x(t)  =  e'^(*“*®^xo  +  /  ^^Bu(r)dr  (A.4) 

Jto 

This  equation  is  the  primary  building  block  for  discretizing  the  system.  Setting  the 

initial  time  equal  to  {k  —  1)T,  where  T  is  the  discretization  step,  yields 

fkT 

Xfc  =  +  /  e'^(*^“’'^Bu(r)dr  (A.5) 
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In  order  to  perform  the  integration  shown  in  Eq.  (A. 5),  total  knowledge  of  u(t)  for 
(k  —  1)T  <  t  <  kT  is  required.  Given  a  sampled  input,  one  has  knowledge  of  u(t)  only  at 
the  endpoints,  so  some  assumption  must  be  made  about  the  form  of  u(t)  in  between. 

To  begin,  let  t  =  kT  —  t,  so  dt  =  -dr,  and  Eq.  (A.5)  becomes 

Xfc  =  e'^^Xjt-i  +  /  e‘^*Bu(fcT  -  t)dt  (A. 6) 

Jo 

A. 2  Assumption  of  Constant  Input  Between  Samples 

In  this  scenario,  it  is  assumed  that  the  input  is  constant  and  equal  to  that  at  the 
initial  time,  j.e., 

u(t)  =  u  ((A:  -  l)r)  =  Uk-i,  {k  —  l)r  <t<kT 

Now,  since  the  input  is  a  constant,  it  can  be  brought  outside  the  integration,  so  Eq.  (A. 6) 
becomes 

e^^dt  Bufc_i  (A. 7) 

Integrating  the  term  in  square  brackets  yields 

e^^dt  =  A-^e^*]  ^  =  A-i  (e^^  -  l)  ( A.8) 

so  provided  that  A  is  invertible, 

Xfc  =  e'^^Xfc_i  +  A“^  (e^^  -  l)  b]  Ufc_i  (A.9) 

Upon  inspection,  it  is  apparent  that  Eq.  (A.9)  has  the  same  form  as  Eq.  (A.2a), 

where 

Ad  =  and  B^  =  A“^  ^e'^^  -  B  (A.  10) 

This  can  be  converted  to  a  discretized  transfer  function  via 

^  =  Cd{zI-Adr^Bd  (A.ll) 

u[z) 
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A. 2.1  Example  System  from  Chapter  V.  The  continuous-time  transfer  function  for 


this  system  is 


T{s)  = 


4.8s  +  1.44 
s2  -1- 0.84s -I-  1.44 


This  can  be  converted  to  the  control  canonical  state  space  form  [28] 


x(t)  = 


0  1  0 

x(f)  -1-  u(t) 

-1.44  -0.84  1 


y{t)  =  1.44  4.8  x(/)  +  0  u(t) 


Using  Eqs.  (A.3)  and  (A. 10),  setting  T  =  0.1,  the  discretized  matrices  are 


Ad  = 


0.9930  0.0957 


-0.1378  0.9126 


0.0049 


0.0957 


Cd  =  1.44  4.8 


Dd  =  0 


Therefore,  using  Eq.  (A. 11),  the  resulting  discretized  transfer  function  is 


T(z)  = 


0.46632  -  0.4525 
22-  1.90562+0.9194 


(A.12a) 


(A.12b) 


(A.13) 


A.3  Assumption  of  Linear  Form  Between  Input  Samples 


In  this  scenario,  it  is  assumed  that  the  input  is  linear  between  the  two  samples. 


Therefore,  the  input  is  given  by 


U(,l)  -  - ^ -  Ufc  +  1 - - -  Ufc-l 

t-(k-l)T, 

=  Ufc-l  +  - ^ — ^(Ufc  -  Ufc_i) 


/ 1  rri  .\  ,  1  J 

u{kT  -t)  =  Uk+ - - - 1 

Substituting  Eq.  (A. 14)  into  Eq.  (A.6)  and  simplifying  yields 


(A.14) 


Xfc  =  e-^^Xfc-i  + 


T  "  1  *  r 

/  e^^dt  Bufc  +  —  /  te^^dt  B(ufe_i  -  u^)  (A.15) 

Jo  J  Jo 


A-3 


Integrating  the  second  bracketed  term  by  parts  yields 


te^^dt  =  A-^te^*]  I-  £  A-^e^Vi  =  -  A'^  (e^'^  -  l)  (A.16) 

Now,  substituting  Eqs.  (A.8)  and  (A.16)  into  Eq.  (A.15), 


where 

Ad  =  (A.17a) 

Bdo  =  (a^t)'"  (e^^  -  l)  -  A-i  B  (A.17b) 

Brfi  =  A-^e^^  -  (a^t)”^  -  l)  B  (A.17c) 

As  one  can  see,  this  system  does  not  have  the  form  of  Eq.  (A.2).  Rather,  it  has  the  form 

Xfe  =  AdXfc_i  +  BdoUfc  +  BdiUfc-i  (A.18a) 

Yk  =  CdXjt_i  +  DdUjt_i  (A.18b) 

This  can  be  rewritten  in  2-domain  form  as 


zx(2)  =  Adx(2)  -1-  2BdoU(2)  +  Bdiu(2) 
y(2)  =  Cdx{z)  +  D</u(2) 


so 

x(2)  =  (2I  -  Ad)"^  (^Bdo  +  Bdi)u(2) 

y(z)  =  Cd  [(2I  -  Ad)”^  (zBdo  +  Bdi)  u(2)]  -|-  Ddw(2) 

=  [Cd  -  Ad)~^  Bdi  -I-  Dd  u(2)  -1-  2  [Cd  (zl  -  Ad)“^  Bdo]  u(2) 

Therefore,  the  transfer  function  for  a  system  described  by  Eq.  (A.  18)  is 

^  =  2  [Cd  (2I  -  Ad)-'  Bdo]  +  [Cd  (zl  -  Ad)-'  Bdi  +  Dd]  (A.19) 
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A. 3.1  Example  System  as  in  Chapter  V.  As  before,  the  continuous-time  state  equa¬ 


tions  are  given  in  Eq.  (A. 12).  The  discrete-time  matrices  in  Eqs.  (A. 3)  and  (A. 17),  for 


T  =  0.1,  are  then  given  by 


Ad  = 


0.9930  0.0957 


-0.1378  0.9126 


0.0016 


0.0486 


Cd  =  1.44  4.8 


B^i  = 


0.0032 


0.0471 


Dd  =  0 


Then,  following  Eq.  (A. 19),  the  transfer  function  for  this  model  is 

,  0.2355^2^0.0028.2-  0.2244 


T(z)  = 


22  -  1.90562-1-0.9194 


(A.20) 


A. 3. 2  Further  Approximations  of  the  Tustin  Transformation.  The  Tustin  transfor¬ 
mation  assumes  the  input  has  the  linear  property  discussed  this  section,  plus  it  uses  the 


following  bilinear  approximation  for  the  exponential: 


1  + 


■)(-f 


Using  this  simplification,  the  discretized  matrices  in  Eqs.  (A.3)  and  (A. 17)  are 


Ad  = 


0.9931  0.0956 


-0.1377  0.9128 


Bdo  =  Bdi  = 


0.0024 


0.0478 


Cd  =  1.44  4.8 


and  the  Tustin  transformed  discrete-time  transfer  function  is  then  given  by 


T(z)  = 


0.23302^  -I-  0.00692  -  0.2261 
22  -  1.90592  4-0.9197 


(A.21) 


(A.22) 


A. 3. 3  Example  from  Chapter  VI.  The  continuous-time  state  space  equations  for 


this  system  are  obtained  from  Eq.  (6.1): 


-0.11  0.1 


0.2  -0.202 


0  0.01 


0.002  0 


1  r  * 

y{t)  =  0  -0.1  x(t)  +  0.1  0  u(0 
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Therefore,  the  appropriate  discretized  matrices  are 


0.9039  0.0865 

r  1 

Ad  = 

Cd  = 

0  ~0.1 

0.1729  0.8244 

0.0432  4.7598 

1 

r 

Bdo  =  Bdi  =  10  ^  * 

11 

Q 

0.1  0 

0.9122  0.4323 

and,  using  Eq.  (A. 19),  the  resulting  transfer  functions  are 

(0.9991^2-  1.7284:2  +  0.7311)10-^ 

= - .2  _  1,7283. +  0.7302  “'W 

(4.3232.22  ^  8.6464^  +  4.3232)10“®  .  . 

^2  -  1.72832  +  0.7302 


or,  using  the  form  of  Eq.  (6.9), 


yi^) 


[(1 


.7293  -  0.7306)22  -  1.72932  +  0.7306]  10 
22  -  1.72832  +  0.7302 

4.3232*10-5(22  +  22  +  1) 
22  -  1.72832  +  0.7302 


— Ul(2) 

[ui(2)  -  U2iz)] 


A. 4  Higher  Order  Assumptions 

The  constant-input  assumption  discussed  in  Section  A. 2  is  perfectly  acceptable  for 
systems  driven  by  computer  controllers,  or  other  systems  where  the  input  to  the  system 
actually  is  held  constant  over  the  sample  time  interval.  This  assumption,  however,  is  not  at 
all  suitable  for  the  identification  of  continuous-time  plants  using  samples  of  the  continuous¬ 
time  input  and  output.  For  that,  one  must  assume  a  higher  order  form  on  the  input. 

The  linear  assumption  discussed  in  Section  A. 3  produces  a  relatively  good  model, 
provided  the  sampling  rate  is  high  enough  to  make  the  input  appear  somewhat  linear 
between  samples.  The  e^^  approximation  made  in  the  Tustin  transformation  does  worsen 
the  system  model  somewhat,  but  it  allows  the  simplification  of  some  common  transfer 
function  forms.  Therefore,  one  can  decrease  the  number  of  needed  parameters,  improving 
the  identification  process. 

If  the  sampling  rate  is  too  slow,  or  the  initial  transient  period  of  the  system  needs  to 
be  used,  then  the  linear/Tustin  assumptions  must  be  abandoned  for  a  higher  order  model. 
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These  higher  order  models  can  assume  a  parabolic  or  sinusoidal  form  for  the  input.  Assum¬ 
ing  these  types  of  inputs,  however,  increases  the  order  of  the  discrete-time  model  above  that 
of  the  continuous-time  model.  This  means  additional,  over-determined,  parameters  need  to 
be  identified.  Care  must  be  exercised  in  this  process,  because  identifying  over-determined 
systems  causes  many  problems. 
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Appendix  B.  Noise  Plots  and  Results  for  Chapter  III 


This  appendix  contains  plots  of  representative  noise,  plots  of  estimated  poles  and 
Bode  plots,  and  numerical  results  for  the  frequency  domain  experiments  in  Chapter  III. 
Figure  B.l  contains  representative  uncorrelated  noise  added  to  the  phasor  quantities.  Fig¬ 
ures  B.4  and  B.5  show  the  representative  uncorrelated  noise  added  to  the  magnitude  and 
phase  angle,  and  the  noise  as  it  appears  on  the  phasor  quantities.  Figures  B.8  and  B.9 
show  the  same  for  the  noise  added  to  the  magnitude  in  dB’s  and  phase  angle  in  degrees. 
The  normalized  numerical  results  for  each  of  the  examined  noise  scenarios  are  given  in  Ta¬ 
bles  B.l,  B.2,  and  B.3.  Shown  are  the  calculated  average  error  in  the  estimate  e  =  0  —  G,  the 
estimation  error  sigma  Oe,  and  the  algorithm  predicted  error  sigma  (Tp  for  the  Monte-Carlo 
analysis.  The  n/a  entries  in  the  tables  indicate  that  these  methods  do  not  produce  an 
estimate  for  the  covariance  value.  The  other  figures  contain  the  final  estimated  Bode  plots 
and  poles  and  zeros  for  the  higher  noise  cases.  Shown  are  the  true  and  average  estimated 
Bode  plots,  along  with  the  one  sigma  bounds  on  the  estimated  Bode  plots.  The  pole  and 
zero  figures  show  each  of  the  estimated  poles  and  zeros,  along  with  the  true  values. 
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Figure  B.l.  Noise  Added  to  A  and  B 
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Table  B.l.  Numerical  Results  for  Noise  on  A  and  B 


Method 

0 

ei%) 

Oe 

(Tp 

e(%) 

Oe 

Cp 

p=0 

p=l 

Least 

Squares 

ai 

-15.0562 

7.7862 

n/a 

-65.0425 

14.9841 

n/a 

ffl2 

-13.8730 

5.7541 

n/a 

-58.5484 

11.4132 

n/a 

Ea 

-14.3308 

6.5080 

n/a 

-60.8256 

12.5748 

n/a 

sa 

-14.1483 

5.9816 

n/a 

-59.3727 

11.5906 

n/a 

Generalized 

Minimum 

Variance 

ffli 

-0.0206 

0.2690 

0.2965 

-0.1449 

0.8501 

0.9362 

0.2 

-0.0265 

0.1816 

0.1764 

-0.1530 

0.5774 

0.5563 

m 

-0.0017 

0.2390 

0.2355 

-0.0520 

0.7538 

0.7441 

62 

-0.0554 

0.5991 

0.6094 

-0.3607 

1.8981 

1.9235 

Generalized 

Least 

Squares 

ffli 

-0.0206 

0.2690 

n/a 

-0.1449 

0.8501 

n/a 

02 

-0.0265 

0.1816 

n/a 

-0.1530 

0.5774 

n/a 

m 

-0.0017 

0.2390 

n/a 

-0.0520 

0.7538 

n/a 

-0.0554 

0.5991 

n/a 

-0.3607 

1.8981 

n/a 

Instrumental 

Variable 

Ol 

0.7340 

5.3380 

11.0606 

3.4853 

17.2366 

11.2673 

02 

-0.1310 

1.4472 

6.1274 

-0.1748 

4.7678 

6.1504 

Ea 

0.0600 

3.2623 

6.6856 

0.7034 

10.4290 

6.7548 

Ea 

■BIMB 

2.3113 

6.2029 

7.5575 

6.2016 

p=2 

p=3 

Least 

Squares 

Ol 

-100.7956 

6.7849 

n/a 

-102.3768 

3.0871 

n/a 

02 

-91.5509 

3.8402 

n/a 

-97.7128 

0.7871 

n/a 

m 

-95.0081 

4.7130 

n/a 

-100.4070 

1.5802 

n/a 

-92.7353 

3.8645 

n/a 

-99.1122 

0.8354 

n/a 

Generalized 

Minimum 

Variance 

Ol 

-1.2510 

2.6522 

2.9257 

-11.6838 

7.1846 

8.3129 

02 

-1.1654 

1.8652 

1.7240 

-9.5793 

6.5824 

4.5607 

Eai 

-0.6375 

2.3487 

2.3345 

-7.3512 

6.8212 

6.8449 

Eai 

-2.9605 

6.0049 

5.9899 

-24.7188 

18.1148 

16.5515 

Generalized 

Least 

Squares 

Ol 

-1.2510 

2.6522 

n/a 

-11.6838 

7.1846 

n/a 

02 

-1.1654 

1.8652 

n/a 

-9.5793 

6.5824 

n/a 

Eai 

2.3487 

n/a 

-7.3512 

6.8212 

n/a 

Eai 

-2.9605 

6.0049 

n/a 

-24.7188 

18.1148 

n/a 

Instrumental 

Variable 

Ol 

-73.9084 

327.3302 

29.3982 

-4224.8129 

40090.8675 

784.6319 

0,2 

-21.2100 

56.3442 

7.3096 

127.2910 

1885.0227 

40.7235 

lai 

-30.0850 

117.6643 

8.7277 

-816.8024 

7149.0848 

144.3396 

m 

-21.3158 

52.7278 

7.3603 

426.4470 

4774.7141 

85.3484 
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Figure  B.4.  Noise  Added  to  Magnitude  and  Radian  Phase 
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Figure  B.5.  Transformed  Noise  On  Magnitude  and  Radian  Phase 


Table  B.2.  Numerical  Results  for  Noise  on  J?  and  (j> 


Method  6 


Least 

Squares 


Generalized  ui 
Minimum  02 
Variance 


Generalized  |  oi 
Least 
Squares 


Variable 


Least 

Squares 


Generalized  ai 
Minimum  02 
Variance 

1  ^2 


Instrumental  a2 
Variable  bi 


e(%) 

(Te 

O'  p 

11 

0 

-13.7859 

8.3089 

n/a 

-12.7156 

6.7777 

n/a 

-13.2254 

7.1137 

n/a 

-13.0119 

7.0774 

n/a 

-0.0333 


-0.0214 


-0.0191 


-0.0253 


-0.0327 


-0.0332 


0.2404 


0.1388 


0.2253 


0.3002 


0.3053 


0.1799 


-0.0810  0.5487 


0.2465 


0.1509 


0.2433 


0.3322 


n/a 


n/a 


n/a 


n/a 


0.5676 

3.6425 

11.0571 

2.3534 

11.3904 

11.2125 

<12 

-0.2263 

1.8593 

6.1239 

-0.2770 

6.4606 

6.1379 

mi 

-0.1548 

1.9778 

6.6750 

0.0437 

6.8177 

6.7071 

mil 

-0.4102 

3.2066 

6.1972 

10.8584 

6.1859 

-99.4138 


-89.8717 


-93.4695 


-91.1127 


-1.0872 


-0.9521 


-0.8351 


-1.4429 


-1.7115 


-1.2550 


-1.0672 


-3.3997 


-46.5619 


-16.0451 


-23.7460 


-16.9293 


p=2 


7.4107 


5.4065 


5.8279 


5.5227 


2.4214 


1.4101 


2.2562 


3.0272 


3.0686 


1.8752 


2.6737 


5.6253 


143.6145 


39.5066 


58.3411 


43.3818 


n/a 


n/a 


n/a 


n/a 


2.4171 


1.4821 


2.4079 


3.2757 


n/a 


n/a 


n/a 


n/a 


20.0687 


6.7231 


7.0748 


6.4407 


e(%) 

o-e 

(Tp 

p=l 

-60.8193 

17.0990 

n/a 

-54.5728 

14.5087 

n/a 

-56.8998 

15.1957 

n/a 

-55.4231 

14.8450 

n/a 

-0.1644 


-0.1232 


-0.1101 


-0.1700 


-0.2092 


-0.1760 


-0.1166 


-0.4561 


0.7582 


0.4388 


0.7114 


0.9492 


0.9694 


0.5730 


0.8457 


1.7458 


0.7779 


0.4763 


0.7685 


1.0490 


n/a 


n/a 


n/a 


n/a 


-102.6223 


-97.3743 


-100.2593 


-98.8350 


-7.6534 


-9.7354 


-7.3839 


-14.5805 


-83.8464 


-78.8845 


-81.7519 


-83.9706 


-187.3172 


-65.2617 


-93.5016 


-63.4651 


p=3 


2.1372 


1.2308 


9.6950 


6.5211 


7.9341 


11.1922 


35.6204 


35.7437 


37.8196 


32.0320 


387.1076 


65.3378 


106.3434 


80.3390 


n/a 


n/a 


n/a 


n/a 


6.5178 


4.0128 


7.0072 


9.0866 


n/a 


n/a 


n/a 


40.0542 


7.0002 


8.7700 


7.7578 
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Frequency  (rad/sec) 


Table  B.3.  Numerical  Results  for  Noise  on  Rib  and  (j)deg 


Method 

9 

e  (%)  (Te 

Cp 

«(%) 

o-e 

Cp 

p=0 

p=l 

Least 

Squares 

ai 

-2.9445 

2.4267 

n/a 

-23.3012 

10.8525 

n/a 

0-2 

-2.8304 

1.6040 

n/a 

-21.0595 

8.6073 

n/a 

m 

-2.9268 

1.7780 

n/a 

-21.9066 

9.1428 

n/a 

-2.9272 

1.8480 

n/a 

-21.4273 

8.8081 

n/a 

Generalized 

Minimum 

Variance 

ffli 

-0.0273 

0.2526 

0.2501 

-0.1481 

0.8007 

0.7902 

0.2 

-0.0190 

0.1429 

0.1546 

-0.1174 

0.4508 

0.4876 

m 

-0.0040 

0.1951 

0.2044 

-0.0010 

0.6165 

0.6463 

-0.0137 

0.2613 

0.3055 

-0.0732 

0.8231 

0.9650 

Generalized 

Least 

Squares 

ai 

-0.0814 

0.4701 

n/a 

-0.4840 

1.5197 

n/a 

«2 

0.2525 

n/a 

-0.1567 

0.8002 

-0.0293 

0.2936 

n/a 

-0.1965 

0.9365 

n/a 

b2 

-0.0719 

0.6230 

n/a 

-0.3082 

1.9785 

n/a 

Instrumental 

Variable 

0.1858 

1.7507 

11.0242 

0.7129 

5.5303 

11.0621 

02 

-0.1124 

0.7662 

6.1226 

-0.2948 

2.4593 

6.1177 

bi 

0.9566 

6.6708 

-0.1318 

3.0493 

6.6748 

b2 

mmmm 

1.2616 

6.2037 

-0.4647 

4.0419 

6.1913 

p=2 

p=3 

Least 

Squares 

ai 

-76.6357 

14.0096 

n/a 

-100.9812 

5.0805 

n/a 

0,2 

-68.3917 

11.0486 

n/a 

-92.4881 

3.6376 

n/a 

la 

-71.2315 

11.7539 

n/a 

-95.8282 

n/a 

la 

-69.2936 

11.1049 

n/a 

-93.6736 

n/a 

Generalized 

Minimum 

Variance 

Oi 

-1.0841 

2.5686 

2.4737 

-9.1244 

8.9164 

7.0716 

02 

-0.9353 

1.4612 

1.5155 

-8.1383 

6.8051 

4.1673 

lai 

0.1100 

1.9513 

2.0448 

6.5124 

6.5088 

lai 

-0.5225 

2.6053 

3.0288 

-4.2140 

9.7562 

9.0579 

Generalized 

Least 

Squares 

o\ 

-3.7523 

5.2088 

n/a 

-40.9184 

28.3573 

n/a 

02 

-0.8867 

2.6295 

n/a 

-14.6615 

28.0427 

n/a 

iai 

-1.6765 

3.1116 

n/a 

-27.5110 

29.5469 

n/a 

b2 

-1.7881 

6.5015 

n/a 

-24.0412 

37.7590 

n/a 

Instrumental 

Variable 

Ox 

1.6252 

28.1418 

11.2921 

21.5836 

1089.3221 

42.1069 

02 

-2.3172 

18.1915 

6.1431 

-29.5097 

81.5201 

7.5469 

bi 

-1.3839 

21.0757 

6.7647 

-12.3400 

236.9992 

10.8943 

wmi 

-2.6517 

21.6585 

6.2209 

-35.4831 

165.4659 

9.4785 
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Appendix  C.  Numerical  Results  for  Chapter  IV 


This  appendix  contains  the  numerical  results  for  the  unmodeled  dynamics  experi¬ 
ments  performed  in  Chapter  IV.  Given  are  the  estimated  parameters  and  the  corresponding 
poles  and  zero  of  the  estimated  plant  for  each  of  the  measurement  windows.  Table  C.l  con¬ 
tains  the  estimated  parameters  for  the  baseline  second-order  plant.  Tables  C.2  and  C.3 
contain  the  estimated  parameters  for  the  fourth-order  plants  containing  the  low  and  high 
frequency  dynamics,  respectively.  Table  C.4  contains  the  estimated  parameters  for  the 
sixth-order  plant.  Finally,  Table  C.5  contains  the  poles  and  zero  for  each  of  the  measure¬ 
ment  windows  in  each  of  the  cases. 
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Table  C.l.  Estimation  Results  Using  Only  the  Second- Order  Plant 


Run 

0-9 

O-e 

(Jp 

9-9 

(Tp 

9 

=  al  =  -0.84 

9 

=  a2  =  -1.44 

1 

1.20  10" 

1.07  10-2 

1.25  10-2 

1.72  10-3 

6.43  10-4 

2 

1.20  10° 

9.82  10-3 

7.95  10-3 

1.44  10° 

4.20  10-3 

1.26  10-3 

3 

1.22  10° 

8.93  10-3 

4.77  10-3 

1.44  10° 

5.34  10-3 

2.28  10-3 

4 

1.20  10° 

1.25  10-1 

9.48  10-3 

1.38  10° 

6.73  10-2 

8.98  10-3 

5 

4.77  10-1 

4.79  10-2 

6.94  10-2 

6.66  10-1 

7.08  10-2 

9.44  10-2 

6 

7.66  10-2 

5.81  10-2 

4.97  10-2 

9.59  10-2 

7.38  10-2 

6.15  10-2 

7 

-1.08  10-3 

9.95  10-3 

1.07  10-2 

-6.00  10-4 

7.46  10-3 

7.01  10-3 

8 

1.08  10-3 

2.37  10-3 

2.70  10-3 

-6.86  10-4 

2.08  10-3 

2.98  10-3 

9 

2.15  10-3 

4.08  10-3 

3.65  10-3 

2.73  10-3 

5.75  10-3 

6.29  10-3 

10 

7.26  10-3 

1.98  10-2 

2.45  10-2 

-7.60  10-4 

1.30  10-2 

1.59  10-2 

11 

4.20  10-2 

1.52  10-1 

1.46  10-1 

3.47  10-2 

8.99  10-2 

1.03  10-1 

12 

-2.95  10-2 

5.88  10-1 

5.19  10-1 

4.05  10-1 

13 

-4.26  10-1 

1.61  10° 

1.86  10° 

1.32  10° 

14 

-7.75  10-1 

3.05  10° 

-9.35  10-3 

3.71  10° 

3.38  10° 

15 

2.39  10-1 

9.67  10° 

5.02  10° 

-3.14  10° 

1.48  lOi 

1.04  lOi 

. . "1 

00 

11 

11 

6 

1  =  62  =  1.44 

1 

-5.16  10° 

9.94  10-3 

—CTfiai 

-1.44  10° 

1.73  10-3 

6.43  10-4 

2 

-5.16  10° 

9.00  10-3 

9.23  10-3 

-1.44  10° 

4.23  10-3 

1.26  10-3 

3 

-5.18  10° 

7.22  10-3 

9.35  10-3 

-1.44  10° 

5.47  10-3 

2.29  10-3 

4 

-5.05  10° 

2.92  10-1 

3.26  10-2 

-1.38  10° 

6.59  10-2 

9.00  10-3 

5 

-2.38  10° 

2.51  10-1 

3.37  10-1 

-6.79  10-1 

7.33  10-2 

9.61  10-2 

6 

-3.71  10-1 

2.87  10-1 

2.38  10-1 

-1.01  10-1 

7.74  10-2 

6.48  10-2 

7 

2.35  10-3 

4.82  10-2 

4.48  10-2 

-9.42  10-4 

1.05  10-2 

9.96  10-3 

8 

-5.17  10-3 

1.29  10-2 

1.70  10-2 

1.18  10-3 

8.79  10-3 

1.11  10-2 

9 

-2.60  10-3 

1.49  10-2 

1.55  10-2 

-1.84  10-2 

3.36  10-2 

3.29  10-2 

-6.96  10-3 

2.24  10-2 

2.02  10-2 

-4.22  10-2 

1.11  10-1 

1.43  10-1 

11 

-6.09  10-^ 

2.28  10-2 

2.34  10-2 

-2.21  10-1 

7.79  10-1 

7.49  10-1 

12 

-4.82  10-3 

1.99  10-2 

2.39  10-2 

1.49  10-1 

3.95  10° 

2.89  10° 

13 

5.72  10-4 

2.43  10-2 

2.37  10-2 

2.08  10° 

1.28  lOi 

7.80  10° 

14 

-5.30  10-3 

2.04  10-2 

2.37  10-2 

3.77  10° 

2.81  lOi 

1.47  lOi 

15 

-1.11  10-2 

3.02  10-2 

2.38  10-2 

-1.04  10° 

4.65  IQi 

2.41  lOi 

Table  C.2.  Estimation  Results  With  the  Low  Frequency  Mode  Addition 


Run 

6-e 

CTe 

(Tp 

e-e 

fe 

e 

=  al  =  -0.84 

e 

=  a2=  -1.44 

1 

8.59  10-1 

2.14  10-4 

1.73  10-4 

1.44  10° 

2.84  10-° 

3.31  10-° 

2 

8.56  10-1 

1.30  10-4 

1.06  10-4 

1.44  10° 

6.08  10-° 

3.22  10-° 

3 

8.38  10-1 

3.75  10-° 

1.91  10-° 

1.43  10° 

6.35  10-4” 

7.31  lO-^" 

4 

8.34  10-1 

3.45  10-° 

2.99  10-° 

1.43  10° 

1.28  10-° 

1.90  10-° 

5 

1.22  10° 

6.46  10-° 

3.12  10-° 

1.45  10° 

1.63  10-° 

5.93  10-4 

6 

3.66  10-1 

4.99  10-° 

4.07  10-° 

3.38  10-4 

7.42  10-° 

5.71  10-° 

7 

-2.75  10-4 

1.01  10-° 

1.08  10-° 

-8.34  10-° 

7.67  10-° 

7.15  10-° 

8 

5.85  10-4 

2.37  10-° 

2.70  10-° 

-1.80  10-4 

2.08  10-° 

2.99  10-° 

9 

2.90  10-° 

4.08  10-° 

3.66  10-° 

5.28  10-° 

5.75  10-° 

6.28  10-° 

10 

9.62  10-° 

1.98  10-° 

2.44  10-° 

1.33  10-° 

1.30  10-° 

1.59  10-° 

11 

4.46  10-° 

1.52  10-4 

1.46  10-4 

3.59  10-° 

9.00  10-° 

1.03  10-4 

12 

-2.75  10-° 

8.01  10-4 

5.88  10-4 

-6.60  10-° 

5.20  10-4 

4.06  10-4 

13 

-4.28  10-1 

2.66  10° 

1.61  10° 

1.78  10-° 

1.86  10° 

1.32  10° 

14 

-7.77  10-1 

5.84  10° 

3.05  10° 

-1.25  10-° 

3.72  10° 

3.39  10° 

15 

2.37  10-1 

9.67  10° 

5.02  10° 

-3.14  10° 

1.48  104 

1.04  104  II 

IHH 

e  =  bl  =  4.8 

e 

=  62  =  1.44 

1 

-4.80  IF^ 

2.71  10-° 

2.47  10-° 

-1.44  10° 

1.29  10-^ 

2 

-4.79  10° 

4.83  10-° 

4.63  10-° 

-1.44  10° 

8.83  10-^ 

3 

-4.79  10° 

4.74  10-° 

6.97  10-° 

-1.44  10° 

1.27  10-° 

5.86  10-° 

4 

-4.80  10° 

8.87  10-° 

1.09  10-4 

-1.44  10° 

2.53  10-° 

1.83  10-° 

5 

-5.48  10° 

5.92  10-° 

4.34  10-° 

-1.46  10° 

2.91  10-° 

1.04  10-° 

6 

-1.35  10° 

2.84  10-4 

2.21  10-4 

-2.52  10-4 

8.67  10-° 

6.62  10-° 

7 

5.21  10-° 

4.96  10-° 

4.60  10-° 

3.75  10-° 

1.10  10-° 

1.03  10-° 

8 

-2.11  10-° 

1.29  10-° 

1.70  10-° 

9.57  10-4 

8.78  10-° 

1.11  10-° 

9 

-1.37  10-° 

1.49  10-° 

1.55  10-° 

-3.46  10-° 

3.35  10-° 

3.28  10-° 

10 

-2.21  10-° 

2.24  10-° 

2.02  10-° 

-6.72  10-° 

1.11  10-4 

1.42  10-4 

11 

-1.63  10-° 

2.27  10-° 

2.34  10-° 

-2.46  10-4 

7.75  10-4 

7.45  10-4 

12 

-2.05  10-° 

1.98  10-° 

2.38  10-° 

1.26  10-4 

3.93  10° 

2.88  10° 

13 

-1.51  10-° 

2.42  10-° 

2.37  10-° 

2.06  10° 

1.28  104 

7.77  10° 

14 

-2.10  10-° 

2.04  10-° 

2.36  10-° 

3.76  10° 

2.80  104 

1.46  104 

15 

-2.68  10-° 

3.01  10-° 

2.38  10-° 

-1.05  10° 

4.64  104 

2.40  104 
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Table  C.3.  Estimation  Results  With  the  High  Frequency  Mode  Addition 


Run 

e-e 

O’e 

(Tp 

9-9 

O-e 

(Tp 

9 

=  al  =  -0.84 

9 

=  a2  =  -1.44 

1 

1.20  10“ 

1.06  10-2 

1.23  10-2 

1.44  10“ 

1.71  10-° 

6.39  10-^ 

2 

1.20  10° 

9.74  10-° 

7.84  10-° 

1.44  10° 

4.18  10-° 

1.25  10-° 

3 

8.96  10-° 

4.72  10-° 

1.44  10° 

5.33  10-° 

2.27  10-° 

4 

1.20  10° 

1.26  10-1 

9.51  10-° 

1.38  10° 

6.80  10-2 

9.01  10-° 

5 

4.74  10-1 

4.82  10-2 

6.97  10-2 

6.66  10-1 

7.10  10-2 

9.44  10-2 

6 

7.28  10-2 

5.82  10-2 

4.98  10-2 

9.69  10-2 

7.36  10-2 

6.14  10-2 

7 

-4.30  10-° 

9.97  10-° 

1.07  10-2 

-5.53  10-° 

7.45  10-° 

7.01  10-° 

8 

1.04  10“° 

2.37  10-° 

2.70  10-° 

-6.10  10-1 

2.08  10-° 

2.98  10-° 

9 

6.13  10-^ 

4.04  10-° 

3.64  10-° 

-9.57  10-° 

5.74  10-° 

6.30  10-° 

10 

-5.64  10-2 

2.16  10-2 

2.54  10-2 

-7.39  10-2 

1.35  10-2 

1.59  10-2 

11 

-1.05  10° 

2.82  10-1 

1.77  10-^ 

-8.38  10-1 

1.14  10-1 

8.86  10-2 

12 

BKinilU 

6.67  10-1 

-7.38  10° 

5.84  10-1 

3.07  10-1 

13 

HiSnO 

5.64  10-2 

-8.95  lOi 

2.84  10-1 

2.91  10-1 

14 

-6.65  10-1 

1.12  10-1 

4.20  10-2 

-2.77  102 

7.10  10-1 

6.96  10-1 

15 

,  .. 

4.03  10-1 

4.47  10-1 

-5.35  102 

5.11  10° 

4.11  10° 

00 

II 

11 

11 

II 

1 

-5.15  10“ 

9.74  10-° 

1.25  10-2 

-1.44  10“ 

1.71  10-° 

6.36  10-1 

2 

-5.16  10° 

8.90  10-° 

9.09  10-° 

-1.44  10° 

4.19  10-° 

1.24  10"° 

3 

-5.17  10° 

7.19  10-° 

9.28  10"° 

-1.44  10° 

5.44  10-° 

2.26  10-° 

4 

-5.04  10° 

2.95  10-1 

3.28  10-2 

-1.38  10° 

6.63  10-2 

8.99  10-° 

5 

-2.37  10° 

2.52  10-1 

3.37  10-1 

-6.82  10-1 

7.31  10-2 

9.56  10-2 

6 

-3.65  10-1 

2.87  10-1 

2.38  10-1 

-1.08  10-1 

7.69  10-2 

6.43  10-2 

7 

9.69  10-° 

4.82  10-2 

4.49  10-2 

-7.16  10-° 

1.05  10-2 

9.91  10-° 

8 

-4.88  10-° 

1.29  10-2 

1.70  10-2 

6.21  10-1 

8.80  10-° 

1.11  10-2 

9 

-8.79  10-° 

1.48  10-2 

1.55  10-2 

4.49  10-2 

3.35  10-2 

3.29  10-2 

10 

-2.23  10-2 

2.22  10-2 

1.97  10-2 

4.58  10-1 

1.23  10-1 

1.49  10-1 

11 

-9.45  10-2 

2.29  10-2 

2.10  10-2 

5.75  10° 

1.43  10° 

8.95  10-1 

12 

-4.53  10-1 

2.41  10-2 

1.79  10-2 

4.47  lOi 

4.67  10° 

3.16  10° 

13 

-3.77  10° 

1.32  10-2 

8.24  10-° 

1.65  lOi 

4.21  10-1 

2.52  10-1 

14 

-4.34  10° 

1.20  10-2 

7.20  10-° 

4.01  10° 

2.82  10-1 

1.67  10-1 

15 

3.06  10-1 

2.68  10-2 

2.04  10-2 

7.32  lOi 

2.44  10° 

2.94  10° 
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Table  C.4.  Estimation  Results  With  Both  Low  and  High  Frequency  Mode  Additions 


Run 

9-9 

Oe 

(Tp 

9-9 

o-e 

dp 

9 

=  al  =  -0.84 

9 

1 

II 

e 

II 

1 

8.59  10-i 

2.14  10-^ 

1.73  lO-'i 

1.44  10*^ 

2.84  10-° 

3.30  10"° 

2 

8.56  10-^ 

1.30  10-^ 

1.06  10-^ 

1.44  10° 

6.08  10-° 

3.22  10-° 

3 

8.38  10-1 

3.75  10-° 

1.91  10-° 

1.43  10° 

6.34  10-° 

7.30  10-° 

4 

8.34  10-1 

3.45  10-° 

2.99  10-° 

1.43  10° 

1.28  10-° 

1.90  10-° 

5 

1.22  10° 

6.41  10-° 

3.09  10-° 

1.45  10° 

1.63  10-° 

5.90  10-“ 

6 

3.60  10-1 

5.03  10-° 

4.11  10-° 

3.36  10-1 

7.43  10-° 

5.73  10“° 

7 

-3.78  10-° 

1.01  10-° 

1.09  10-° 

-7.98  10-° 

7.66  10-° 

7.15  10-° 

8 

5.47  10-^ 

2.37  10-° 

2.70  10-° 

-1.86  10-“ 

2.08  10-° 

2.99  10-° 

9 

1.36  10-° 

4.05  10-° 

3.64  10-° 

-7.12  10-° 

5.74  10-° 

6.29  10-° 

10 

-5.39  10-° 

2.15  10-° 

2.53  10-° 

-7.18  10-° 

1.35  10-° 

1.59  10-° 

11 

-1.04  10° 

2.81  10-1 

1.76  10-1 

-8.38  10-1 

1.14  10-1 

8.86  10-° 

12 

-9.08  10° 

9.90  10-1 

6.66  10-1 

-7.39  10° 

5.86  10-1 

3.07  10-1 

13 

-5.74  10° 

1.37  10-1 

5.64  10-° 

-8.96  lOi 

2.84  10-1 

2.91  10-1 

14 

-6.66  10-1 

1.12  10-1 

4.20  10-° 

-2.77  10° 

7.10  10-1 

6.96  10-1 

15 

-1.44  lOi 

4.03  10-1 

4.47  10-1 

-5.35  10° 

5.11  10° 

4.11  10° 

9  =  bl  =  4.8 

9 

=  62  =  1.44 

1 

-4.80  10° 

2.70  10-° 

2.45  10-° 

-1.44  10° 

1.15  10"° 

1.28  10-° 

2 

-4.79  10° 

4.80  10-° 

4.60  10-° 

-1.44  10° 

1.83  10-° 

8.78  10-° 

3 

-4.79  10° 

4.71  10-° 

6.93  10-° 

-1.44  10° 

1.26  10-° 

5.83  10-° 

4 

-4.80  10° 

8.86  10-° 

1.08  10-“ 

-1.44  10° 

2.52  10-° 

1.82  10-° 

5 

-5.48  10° 

5.81  10-° 

4.28  10-° 

-1.46  10° 

2.90  10-° 

1.03  10-° 

6 

-1.33  10° 

2.85  10-1 

2.22  10-1 

-2.55  10-1 

8.64  10-° 

6.60  10-° 

7 

6.07  10-° 

4.96  10-° 

4.60  10-° 

3.10  10-° 

1.10  10"° 

1.03  10-° 

8 

-1.79  10-° 

1.29  10-° 

1.70  10-2 

2.25  10-“ 

8.79  10-° 

1.11  10-° 

9 

-1.99  10-° 

1.47  10-° 

1.55  10-° 

2.85  10-° 

3.34  10-° 

3.28  10-° 

■■ 

-3.74  10-° 

2.22  10-° 

1.97  10-° 

4.30  10-1 

1.23  10-1 

1.48  10-1 

11 

-1.10  10-1 

2.28  10-° 

2.10  10-° 

5.70  10° 

1.42  10° 

8.90  10-1 

12 

-4.68  10-1 

2.40  10-° 

1.79  10-° 

4.45  lOi 

4.66  10° 

3.15  10° 

13 

-3.77  10° 

1.31  10-° 

8.21  10-° 

1.64  lOi 

4.20  10-1 

2.51  10-1 

14 

-4.34  10° 

1.19  10-° 

7.17  10-° 

4.00  10° 

2.81  10-1 

1.67  10-1 

15 

2.90  10-1 

2.68  10-° 

2.03  10-° 

7.30  lOi 

2.43  10° 

2.93  10° 

Table  C.5.  Estimated  Poles  and  Zeros  for  Each  Case 


Run 

2nd  Order 

4th  Order  Low 

4th  Order  High 

6th  Order 

Poles  -  True  =  — 

0.42-|-il.l24 

1 

0.356  &  0.002 

0.009+j0.032 

0.353  k  0.002 

0.009+j0.032 

2 

0.356  &  0.004 

0.008-fj0.053 

0.354  k  0.004 

0.008-Hj0.053 

3 

0.365  k  0.012 

-0.001-|-j0.074 

0.362  k  0.012 

-0.001-|-j0.074 

4 

0.181+j0.152 

-0.003-!-j0.073 

0.179+j0.157 

-0.003-|-j0.073 

5 

-0.181+j0.861 

0.414  &  -0.031 

-0.183fj0.861 

0.412  k  -0.030 

6 

-0.382+jl.095 

-0.237-l-jl.022 

-0.384+j  1.094 

-0.240-|-jl.023 

7 

-0.421+jl.l24 

-0.420-l-jl.128 

-0.422-hjl.l23 

-0.422-hjl.l27 

8 

-0.419+jl.l25 

-0.420-hjl.l24 

-0.419-|-jl.l25 

-0.420+jl.l24 

9 

-0.419+jl.l23 

-0.419-|-jl.l22 

-0.420+jl.l28 

-0.419fjl.l28 

10 

-0.416+jl.l26 

-0.415-t-jl.l25 

-0.448-fj  1.146 

-0.447+j  1.145 

11 

-0.399+jl.ll6 

-0.398+jl.ll6 

-0.943-hjl.l78 

-0.941+jl.l80 

12 

-0.435+jl.l48 

-0.434-l-jl.148 

-8.934  k  -0.987 

-8.928  k  -0.989 

13 

-0.633+jl.017 

-0.634-hjl.018 

-3.288-fj8.954 

-3.288+j8.955 

14 

-0.808+j0.893 

-0.809-fj0.894 

-0.753+jl6.669 

-0.753-t-jl6.669 

15 

-0.301+j2.119 

-0.301-l-j2.119 

-7.622-t-j21.881 

-7.622+j21.881 

Zero  -  True 

=  0.300 

1 

0.006 

0.006 

2 

0.024 

0.024 

3 

0.149 

0.149 

4 

5.038 

5.596 

5 

-0.035 

-0.033 

6 

-0.302 

-0.344 

-0.300 

-0.341 

7 

-0.300 

-0.305 

-0.298 

-0.303 

8 

-0.301 

-0.300 

-0.300 

-0.300 

9 

-0.296 

-0.294 

-0.310 

-0.307 

10 

-0.292 

-0.287 

-0.397 

-0.393 

11 

-0.254 

-0.250 

-1.528 

-1.521 

12 

-0.331 

-0.328 

-10.602 

13 

-0.732 

-0.732 

-17.372 

-17.372 

14 

-1.087 

-1.087 

-11.838 

-11.838 

15 

-0.082 

-0.082 

-14.625 

-14.624 

Appendix  D.  Results  for  Chapter  V 

This  appendix  contains  the  normalized  numerical  results  for  the  higher  noise  exper¬ 
iments  in  Chapter  V.  Tables  D.l  and  D.2  contain  the  results  for  experiments  with  input 
noise  only,  and  the  others  contain  the  results  for  the  experiments  with  both  input  and  out¬ 
put  noise.  Given  in  each  case  are  the  actual  average  error  (e)  and  error  sigma  (<Te),  along 
with  the  algorithm  predicted  error  sigma  {(Tp)  for  the  100  run  Monte-Carlo  analysis.  The 
NaN  and  Inf  entries  indicate  that  the  particular  algorithm  did  not  converge  in  that  case. 
NaN  and  Inf  are  the  not-a-number  and  infinity  values,  respectively,  that  are  produced  by 
Matlab. 


Table  D.l.  Numerical  Results  for  p  =  3  -  Output  Noise  Only 


Method 

e 

e  =  e~6{%) 

o-e 

(Jp 

9-e{%) 

Oe 

(7p 

Initial  Time 

No  Initial  Time 

Least 

Squares 

ai 

-16.3283 

1.3336 

0.8827 

-15.4486 

1.2063 

0.8573 

02 

-30.5571 

2.5130 

1.6678 

-29.0724 

2.2751 

1.6293 

El 

-11.8435 

1.2978 

6.1293 

-9.0016 

1.2583 

6.2259 

-34.4039 

2.7767 

6.9777 

-29.2986 

2.6407 

6.9806 

Generalized 

Minimum 

Variance 

Oi 

0.0005 

0.0381 

0.0413 

-0.0034 

0.0404 

0.0415 

0.2 

0.0016 

0.0825 

0.0882 

-0.0066 

0.0871 

0.0885 

h 

-0.0244 

0.3395 

0.3443 

0.0309 

0.3428 

0.3444 

Ell 

-0.0268 

0.3392 

0.3470 

0.0301 

0.3391 

0.3467 

Generalized 

Least 

Squares 

Ol 

0.0002 

0.0380 

0.0404 

-0.0022 

0.0407 

0.0417 

02 

0.0011 

0.0819 

0.0862 

-0.0040 

0.0876 

0.0891 

Ell 

-0.0227 

0.3358 

0.3431 

0.0265 

0.3420 

0.3422 

Ell 

-0.0251 

0.3374 

0.3461 

0.0267 

0.3378 

0.3444 

Instrumental 

Variable 

Ol 

-0.0003 

0.0422 

0.9753 

0.0038 

0.0454 

0.9410 

02 

-0.0011 

0.0855 

1.8659 

-0.0038 

0.0895 

1.8151 

Ell 

-0.0128 

1.1521 

6.1058 

0.2370 

0.8167 

6.2204 

m 

-0.0038 

1.1974 

6.8403 

0.3110 

0.9930 

6.8564 
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Table  D.2.  Numerical  Results  for  p  =  4  -  Output  Noise  Only 


Method 

Oe 

C  p 

Oe 

(Tp 

Initial  Time 

No  Initial  Time 

ai 

-80.2695 

1.8686 

2.0634 

-79.8178 

1.7720 

2.0549 

Least 

0.2 

-133.0600 

3.9315 

4.0202 

-132.8790 

3.7436 

4.0182 

Squares 

m 

-217.6903 

32.6668 

59.1150 

-213.0996 

35.6261 

60.1482 

-362.5656 

37.2700 

64.9588 

-350.6757 

40.2495 

65.5272 

Generalized 

Ol 

-0.0175 

0.4260 

0.3848 

-0.0771 

0.4257 

0.3857 

Minimum 

02 

-0.0171 

0.9154 

0.8188 

-0.1460 

0.9103 

0.8215 

Variance 

m 

0.1072 

3.3889 

3.3907 

0.6370 

3.4707 

3.3981 

WM 

0.0526 

3.3985 

3.4191 

0.6210 

3.4506 

3.4228 

Generalized 

Ol 

0.0138 

0.4245 

0.3764 

-0.0271 

0.4309 

0.3850 

Least 

02 

0.0487 

0.9086 

0.8003 

-0.0425 

0.9199 

0.8199 

Squares 

m 

0.0009 

3.3775 

3.3745 

0.4608 

3.4922 

3.3691 

m 

-0.0650 

3.4030 

3.4051 

0.4416 

3.4689 

3.3922 

Ol 

NaN 

NaN 

NaN 

0.0259 

0.4596 

9.1528 

Instrumental 

02 

NaN 

NaN 

NaN 

-0.0614 

0.9135 

17.7292 

Variable 

lai 

NaN 

NaN 

NaN 

2.4422 

8.1945 

61.4004 

62 

NaN 

NaN 

NaN 

3.1748 

9.9608 

67.2974 

Table  D.3.  Numerical  Results  for  p  =  3  -  Input  and  Output  Noise 


Method 

e 

Oe 

CTp 

9-0{%) 

(Xp 

Initial  Time 

No  Initial  Time 

Least 

Squares 

Ol 

-17.2891 

1.2256 

0.8709 

-15.9980 

1.1173 

0.8495 

02 

-31.2484 

2.3315 

1.6623 

-29.0062 

2.1094 

1.6277 

■a 

-43.0761 

7.7805 

4.7164 

-39.4903 

7.5683 

4.7598 

WM 

-69.7460 

8.8178 

5.3556 

-62.9591 

8.5739 

5.3315 

Generalized 

Minimum 

Variance 

Ol 

0.0015 

0.0877 

-0.0130 

0.0899 

0.0873 

02 

-0.0053 

bi 

-0.5206 

0.6930 

0.7201 

-0.4320 

0.6973 

0.7155 

mm 

-0.5089 

0.7059 

-0.4174 

0.7041 

0.7237 

Generalized 

Least 

Squares 

Ol 

0.0328 

0.0404 

IMiIilW 

0.0929 

0.0417 

02 

0.0576 

0.2029 

0.0860 

0.1965 

0.0891 

hi 

-0.7842 

0.7207 

0.3419 

-0.7059 

0.7527 

0.3413 

lai 

-0.7752 

0.7381 

0.3449 

-0.6906 

0.7611 

0.3434 

Instrumental 

Variable 

Ol 

-7.0318 

2.6997 

1.2022 

-5.4413 

2.3233 

1.1499 

-8.7355 

3.5016 

2.0832 

-6.3069 

2.8940 

2.0245 

lai 

-60.5149 

17.1172 

5.5997 

-55.0791 

16.3389 

5.5930 

m 

23.9440 

7.0957 

-69.8840 

22.1776 

6.8645 
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Table  D.4.  Numerical  Results  for  p  =  4  -  Input  and  Output  Noise 


Method 

9 

IKH 

(T  p 

lEB^I 

(Tp 

Initial  Time 

No  Initial  Time 

Least 

Squares 

ai 

-76.0406 

2.0994 

1.9559 

-75.8311 

1.8472 

1.9655 

a2 

-135.6182 

4.3837 

3.9279 

-135.1656 

3.9999 

3.9403 

wm 

-65.5936 

10.0666 

8.8508 

-63.9891 

10.2143 

8.7801 

-142.2788 

10.7713 

9.3633 

-138.8498 

10.6890 

9.2873 

Generalized 

Minimum 

Variance 

ai 

-0.1503 

0.5865 

0.6404 

-0.3038 

0.6772 

0.6487 

02 

-0.7400 

1.2970 

1.3586 

-1.0400 

1.4508 

1.3771 

h 

-28.3322 

4.2703 

4.8447 

-27.2686 

4.7292 

4.8908 

wm 

-27.5265 

4.2037 

4.9049 

-26.4482 

4.7150 

4.9586 

Generalized 

Least 

Squares 

ai 

NaN 

NaN 

Inf 

NaN 

NaN 

Inf 

02 

NaN 

NaN 

Inf 

NaN 

NaN 

Inf 

iai 

NaN 

NaN 

Inf 

NaN 

NaN 

Inf 

m 

NaN 

NaN 

Inf 

NaN 

NaN 

Inf 

Instrumental 

Variable 

Ol 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

02 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 

IBI 

NaN 

NaN 

NaN 

NaN 

NaN 

m 

NaN 

NaN 

NaN 

NaN 

NaN 

NaN 
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